NAT'L  INST.  OF  STAND  & TECH  R.I.C 


Penetration  of  Proton  Beams 
Through  Water 

II.  Three-Dimensional  Absorbed 
Dose  Distributions 


I 


Martin  J.  Berger 


U.S.  DEPARTMENT  OF  COMMERCE  ^ 

Technology  Administration  , 

National  Institute  of  Standards 

' ' ’C  ' ' 

and  Technology  -h 

Gaithersburg,  MD  20899 

Prepared  fon  p 

National  Cancer  Institute 

Bethesda,  MD  20892  ffiS . m 


■■ 

. 


i 


I 


1 

\ 


i 

1 

i 


Penetration  of  Proton  Beams 
Through  Water 

II.  Three-Dimensional  Absorbed 
Dose  Distributions 


Martin  J.  Berger 


U.S.  DEPARTMENT  OF  COMMERCE 
Technology  Administration 
National  Institute  of  Standards 
and  Technology 
Gaithersburg,  MD  20899 

Prepared  for; 

National  Cancer  Institute 
Bethesda,  MD  20892 

December  1993 


U.S.  DEPARTMENT  OF  COMMERCE 
Ronald  H.  Brown,  Secretary 

TECHNOLOGY  ADMINISTRATION 

Mary  L Good,  Under  Secretary  for  Technology 


NATIONAL  INSTITUTE  OF  STANDARDS 
AND  TECHNOLOGY 
Arati  Prabhakar,  Director 


Penetration  of  Proton  Beams  Through  Water 
n.  Three-Dimensional  Absorbed  Dose  Distributions 


Martin  J.  Berger* 


Physics  Laboratory** 

National  Institute  of  Standards  and  Technology 
Gaithersburg,  MD  20899 


Abstract 

This  report  describes  methods  and  computer  programs  for  calculating 
absorbed-dose  distributions  in  a water  target  irradiated  by  proton  beams.  The 
spatial  pattern  of  absorbed  dose  from  monoenergetic  pencil  beams  is  calculated  as 
a function  of  the  depth  in  the  target  and  of  the  radial  distance  from  the  pencil 
beam.  This  calculation  uses  the  Monte  Carlo  program  PTRAN  and  also  Moli^re’s 
theory  of  radial  multiple-scattering  deflections.  Such  pencil-beam  results  are  then 
combined  linearly  to  obtain  three-dimensional  absorbed-dose  distributions  from 
beams  that  irradiate  a field  with  arbitrary  sh^e  and  size.  In  general,  this  requires 
a double  numerical  quadrature.  Specialized  formulas  for  circular  and  rectangular 
fields  are  also  given  which  require  only  a single  numerical  quadrature. 


*Work;  carried  out  for  the  National  Institute  of  Standards  and  Technology  under  contract  50SBNB3C7583. 

**Ionizing  Radiation  Division,  National  Institute  of  Standards  and  Technology,  Technology 
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1.  Introduction 


This  is  the  third  report  in  a series  dealing  with  the  transport  of  proton  beams  through 
water.  The  first  report  (Berger,  1993a)  described  the  Monte  Carlo  transport  program  PTRAN 
which  calculates  the  penetration,  diffusion  and  slowing  down  of  protons  in  an  extended  medium. 
The  second  report  (Berger,  1993b)  presented  depth-dose  curves,  proton  spectra  and  LET 
distributions  as  functions  of  the  depth  in  a water  medium. 

The  present  report  deals  with  the  calculation  of  absorbed  dose  distributions  as  a function 
of  two  or  three  spatial  variables.  Section  2 discusses  methods  for  obtaining  dose  distributions 
from  beams  with  arbitrary  cross  sections  from  the  superposition  of  dose  distributions  from  pencil 
beams.  In  general  this  requires  a double  numerical  integration.  However,  for  beams  with 
circular  or  rectangular  cross  sections  one  of  the  integrations  is  done  analytically,  so  that  only  a 
single  numerical  quadrature  is  required.  In  section  3 a database  of  depth-dose  distributions  is 
described  which  was  obtained  widi  the  Monte  Carlo  program  PTRAN,  and  which  gives  the 
depth-dependence  for  beams  with  initial  energies  from  250  MeV  to  50  MeV.  In  section  4 the 
radial  dependence  of  the  absorbed  dose  from  pencil  beams  is  discussed.  It  is  shown  that  the 
radial  dependence  obtained  with  the  PTRAN  program  is  in  good  agreement  with  the  dependence 
that  can  be  obtained  more  directly  and  simply  with  Molibre’s  theory  of  radial  multiple-scattering 
deflections.  However,  at  depths  equal  to  or  greater  than  the  depth  where  the  Bragg  peak  occurs, 
it  is  still  necessary  to  rely  on  the  Monte  Carlo  method.  In  section  5,  various  examples  are  given 
of  absorbed-dose  distributions  from  beams  with  circular  and  rectangular  cross  sections.  In 
section  6,  computer  programs  are  described  with  which  the  database  of  depth-dose  distributions 
can  be  accessed,  and  with  which  the  superposition  of  pencil-beam  results  can  be  accomplished. 


2.  Superposition  of  Pencil  Beam  Results 

We  consider  a narrow  pencil  beam  incident  along  the  z axis  onto  a water  medium  that 
occupies  the  region  z > 0.  It  is  convenient  to  represent  the  absorbed-dose  distribution  from  such 
a beam  in  terms  of  two  quantities:  the  energy-deposition  distribution  dD/dz,  and  the  radial  dose 
distribution  f(p,z),  at  depth  z.  These  distributions  are  defined  as  follows: 

dz  is  the  amount  of  energy,  per  incident  proton,  which  is  imparted  to  the 

^ medium  at  depths  between  z and  z + dz. 

2x-p  f(p,z)dp  is  the  fraction  of  the  energy  (dD/dz)dz  that  is  imparted  to  the  medium  at 
radial  distances  between  p and  p + dp  from  the  z-axis. 

In  the  present  work  dD/dz  was  obtained  with  the  Monte  Carlo  Program  PTRAN, 
combined  with  an  estimate  of  energy  deposition  by  secondary  charged  particles  from  nuclear 
reactions  (Berger,  1993b).  The  calculation  of  dD/dz  will  be  discussed  further  in  section  3,  and 
the  calculation  of  f(p,z)  in  section  4. 

We  now  consider  an  incident  proton  beam  that  consists  of  a bundle  of  parallel  pencil 
beams  of  the  same  energy  that  irradiate  a specified  field  on  the  surface  of  the  phantom.  Within 
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this  field,  let  n(x,y)dxdy  be  the  number  of  pencil  beams  incident  between  x and  x + dx  and 
between  y and  y + dy.  Let  Dg(x,y,z)  denote  the  absorbed  dose  at  a point  (x,y,z),  where  z 
indicates  the  depth  and  x and  y indicate  the  lateral  coordinates.  D^(x,y,z)  can  be  calculated  as  a 
linear  combination  of  dose  distributions  from  pencil  beams: 

Da(x,y,z)  = [ ‘**'1  , (1) 

where 

p = [(x-x')2  + ® 


The  double  integral  in  eq  (1)  must  be  extended  over  the  field  irradiated  by  the  proton 
beam.  With  the  density  n(x,y)  in  units  of  cm“^,  and  dD/dz  in  MeV  cm^/g,  has  units  of 
MeV/g.  (Note  that  1 MeV/g  corresponds  to  1.6022  x 10”^®  Gy).  A computer  program  FREC2 
is  described  in  section  6.6  in  which  eq  (1)  is  applied  to  a rectangular  field. 

It  is  often  advantageous  to  switch  from  Cartesian  to  cylindrical  coordinates,  especially 
when  one  is  considering  a constant  beam  density  n(x,y)  = n^.  Depending  on  the  shape  of  the 
irradiated  field,  it  may  be  then  be  possible  to  integrate  analytically  with  respect  to  the  angular 
variable,  so  that  only  a single  numerical  quadrature  in  the  radial  variable  is  required. 


Figure  1 shows  a curve  in  the  plane  z = 0 which  represents  the  boundary  of  the 
irradiated  field.  Also  shown  is  a circle  of  radius  p around  the  point  P = (x,y,0).  Let  yl/Jp)  be 
the  fraction  of  the  circular  arc  that  lies  within  the  irradiated  field.  Assuming  a constant  density 
n , eq  (1)  can  be  restated  as 


Da  = no  ^ F^(Z) 


(3) 


where 


FredW  = 2x  f 

J 


00 


f(p,z)  ^p(p)  p dp 


(4) 


For  a uniformly  irradiated  field  F^^j,  is  always  smaller  than  unity.  We  shall  call  it 
reduction  factor^  because  it  represents  the  reduction  of  the  absorbed  dose  compared  to  that  which 
would  prevail  if  the  irradiated  field  were  unbounded.  The  reduction  factor  can  evaluated 
analytically  for  simple  fields  sh^es  such  as  a circular  or  rectangular  field. 


2.1  Circular  Fields. 

Assume  that  the  irradiated  field  is  bounded  by  a circle  with  radius  R,  and  consider  a point 
P = (x,y,0)  located  at  a distance  r from  the  center  of  the  field.  The  reduction  factor  at  point  P 
can  be  shown  to  be  given  by  the  following  formulas: 


2 


For  r < R point  (P  inside  the  field), 


Fred(z.R.r)  = 2x 
and  for  r > R 

Fred(z,R,r)  = 


f R-r  f R+r 

I f(p,z)  p dp  + 27r  I f(p,z)  i/'(p,R,r)  p dp  , 
Jo  J R-r 


f r+R 

[ f(p,z)  \('(p,R,r)  p dp 
Jr-R 


(5a) 

(5b) 


where 


0(p,R,r)  = — cos 

TT 


-1 


2p  r 


1/2 


(6) 


A computer  program  FCIR  is  described  in  section  6.5,  which  calculates  reduction  factors 
and  absorbed-dose  distributions  with  eqs  (5)  and  (6). 

If  the  field  is  not  uniform  but  can  be  represented  by  a density  n(r'),  where  r'  is  the 
distance  from  the  center  of  the  circular  field,  the  absorbed  dose  must  be  calculated  a double 
integral.  When  r < R, 


F^(z,R,r)  = f f(p,z)  p dp  f n[r'((i>)]  d<p  + f f(p,z)  p dp  f n[r'(v>)]  d-p  , (7a) 
Jo  J -X  J R-r  J 


and  when  r > R 


Fred(z>I^>'') 

with 

and 


f f(p,z)  p dp  f n[r'((»)]  d<p 
J r-R  J -<Pi 


■’  = (r^  - 2 


rp  cos^  + p 


2 15  2 


2)1/2  ^ 


<Pl  = COS 


-1  2 + r - R 
2p  r 


(7b) 


(8) 


(9) 


2.2  Rectangular  Fields 

As  was  shown  by  Meredith  and  Neary  (1944),  a point  inside  or  outside  a rectangular  field 
can  be  considered  as  being  at  the  comer  of  four  equivalent  rectangular  fields  which,  when  added 
and/or  subtracted,  represent  the  original  field  (see  fig.  2).  As  a preliminary,  we  therefore 
consider  a circle  with  radius  p centered  on  a comer  of  a rectangle  with  sides  s^  and  S2.  Let 
a = max  (Si,S2)  and  b = min(Si,S2),  and  let  i/'^.(p,Si,S2)  denote  the  fraction  of  the  circle  that  lies 
within  the  rectangle.  It  can  be  shown  that  this  fraction  is  given  by  the  following  expressions: 
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I 


^c(P»Sl,S2)  = ~ 


1 -1 

4 

P 

1 -1 

b 

-1 

a 

= — - cos 

“COS 

4 

P 

P 

= 0 


if  p < b , 
if  b < p < a , 

(10) 

if  a < p <(a^  + , 

if  p > (a2  + 


Next  we  consider  a uniformly  iiradiated  rectangular  field  widi  sides  2 A and  2B, 
occupying  die  region  -A  < x < A and  -B  < y < B.  Hie  fraction  of  a circle  of  radius 
around  the  point  P = (x,y,0)  that  lies  within  ti^e  rectangular  field  can  be  obtained  by  adding 
and/or  subtracting  the  values  of  at  the  comers  of  four  equivalent  rectangles.  The  four  comer 
functions  needed  are 

'Pci  = if'cO».A+x,B+y)  , 

'Pc3  = lAc(P.A+x,|B-y|)  , 


'Pc2  = 4'c(PA-^-x\yB*y)  . 

'P(A  = 'Pc(f>’\^-^\’\^-y\')  ■ 


Four  cases  must  be  considered,  depending  on  die  location  of  P: 

1)  If  X < A and  y < B (P  inside  the  rectaigle),  dien  the  reduction  factor  is 

F^(z,A,B,x,y)  = 2x  [ f(p,z)  p dp 

J O 

+ 2x  [ ^ f(p,z)  + iAc2  + 'Pc3  ^ 'Poi)  p ^ . 

Jpi 

where  p^  = min(A-x,B-y)  andp2  — [(A+x)^  + (B+y)^]^^^. 

2)  If  X < A and  y > B,  then 

Fred(z.A,B,x,y)  = 2x  f ^ f(p,z)  ” Pc3  “ 'P<a)  P . 

where  p^  = y—B  and  P2  = [(A+x)^  + (B+y)^]^^^. 

3)  If  X > A and  y < B,  then 

F^(z,A,B,x,y)  = 2ir  f(p,z)  - 4^2  * Pc3  - 'Pca)  P ^ > 

Jpi 

where  pj  = x— A and  P2  = [(A+x)^  + (B+y)^]^^. 


(12a) 


(12b) 


(12c) 
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(12d) 


4)  If  X a A and  y a B,  then 

Fred(z.A,B,x,y)  = 2i  [ t{p,z)  (i^d  - 'P^2  - + '('c4)  P dp  , 

> Pi 

where  = [(x— A)^  + (y— and  P2  = [(A+x)^  + (B+y)^]^^^. 

A computer  program  FREC  is  described  in  section  6.6  that  calculates  reduction  factors 
and  absorbed-dose  distributions  for  rectangular  fields  using  eqs  (10)  to  (12). 


3.  Database  of  Depth-Dose  Distributions 


The  Monte  Carlo  program  PTRAN  (Berger,  1993a)  was  used  to  calculate  the  following 
two  quantities  that  describe  the  rate  of  energy  loss  of  the  primary  protons  per  unit  depth:  the  loss 
due  to  Coulomb  interactions,  (dE/dz)^.,  and  die  energy  loss  due  to  nuclear  reactions,  (dE/dz)^. 
These  were  evaluated  for  monoenergetic  beams  incident  with  25  energies  between  250  MeV  and 
50  MeV,  at  a grid  of  43  scaled  depths  z/r^  between  0 and  1.04,  where  z is  the  actual  depth  and 
r^  the  CSDA  range.  By  expressing  all  of  the  quantities  in  eq  (13)  as  functions  of  z/r^  rather  than 
r^,  they  become  slowly  varying  functions  of  T^,  which  in  facilitates  interpolation  with  respect  to 
Tq.  Table  1 gives  a set  of  pertinent  range  values  which  were  calculated  with  program  PSTAR 
(Berger,  1992). 


Following  the  procedure  discussed  in  Berger  (1993b)  the  energy  deposited  in  the  phantom 
per  unit  depth  was  calculated  from  the  expression 


dz 


dz 


+ a(z/r„,T„) 


(13) 


where  T^  is  the  beam  energy  and  is  an  absorption  factor  that  takes  into  account  the 

energy  imparted  by  to  the  medium  by  secondary  charged  particles  from  nuclear  reactions.  The 
values  of  lie  between  0.7  and  0.3  for  T^  between  250  MeV  and  50  MeV  and  for  z/r^ 

between  0 and  l.()4.  Values  of  dizh^^T^  were  derived  from  a combination  of  PTRAN  transport 
results  with  of  calculations  of  nonelastic  nuclear  interactions  by  Seltzer(1993).  Tables  of 
a(z/rQ,TQ)  can  be  found  in  Berger  (1993b). 


For  each  energy  T^,  a set  of  1 million  Monte  Carlo  histories  of  primary  protons  was 
sampled  and  analyzed.  Each  set  of  histories  was  divided  into  10  groups  of  100,000  histories,  to 
provide  a basis  for  estimating  the  statistical  error.  The  relative  standard  deviation  of  dD/dz  (ratio 
of  the  standard  deviation  to  die  mean  value)  at  a given  value  of  z/r^  was  found  to  be  nearly 
independent  of  T^.  The  relative  standard  deviations  as  functions  of  the  depth  are  estimated  to 
have  the  following  values: 
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Relative 

z/r^  Standard  Deviation 

of  dD/dz 


< 0.9 

0.01  % 

0.95 

0.02  % 

1.0 

0.1  % 

1.01 

0.3  % 

1.02 

1 % 

1.03 

3 to  6% 

Except  at  depths  comparable  with  and  beyond  the  CSDA  range,  the  statistical  errors  are 
small  compared  to  the  systematic  errors  of  dD/dz.  The  systematic  errors,  estimated  to  be  2 to 
3 percent,  are  due  to  the  uncertainties  of  the  stopping  powers  and  nuclear  interactions  cross 
sections,  and  due  to  approximations  made  in  the  PTRAN  Monte  Carlo  model. 

Curves  of  dD/dz  z/r^  ^e  shown  in  figure  3 for  seven  beam  energies  between  250 
MeV  and  50  MeV.  It  can  be  seen  that  these  curves  have  rather  similar  shapes,  so  that  accurate 
interpolation  with  respect  to  is  possible.  Two  computer  programs  were  developed  which 
carry  out  such  an  interpolation,  using  as  a database  a set  of  curves  of  dD/dz  vs  z/r^  for  25  values 
of  Tq.  Program  PTPOL  calculates  dD/dz  for  a monoenergetic  beam  of  specified  energy. 

Program  PTGPOL  calculates  dD/dz  for  a beam  with  a Gaussian  energy  spectrum  with  specified 
mean  energy  and  standard  deviation.  Descriptions  of  these  programs  are  given  in  the  Sections 

6.1  and  6.2. 


4.  Radial  Dose  Distributions  from  Pencil  Beams 

4.1  Monte  Carlo  Results 

Radial  distributions  from  PTRAN,  for  T^  = 160  MeV  and  70  MeV,  are  shown  in 
figure  4.  Each  of  this  distribution  is  based  on  the  results  from  a sample  of  10  million  proton 
Monte  Carlo  histories.  It  is  useful  to  scale  the  radial  distributions,  by  plotting  the  dimensionless 
quantity  Ittp  z f(p,z)  as  a function  of  the  ratio  p/z.  Figure  4a  shows  scaled  radial  distributions  at 
depths  z = 0.1,  0.3,  0.5  and  0.7  r^,  and  figure  4b  at  z = 0.99,  1,  1.01  and  1.02  r^,  at  and 
beyond  the  Bragg  peak.  As  a result  of  scaling,  the  distributions  at  160  MeV  (solid  curves)  are 
very  similar  to  those  at  70  MeV  (dotted  curves). 

4.2  Results  From  Other  Methods 

A simple  and  convenient  theory  of  the  radial  distribution  of  absorbed  dose  from  proton 
beams  was  first  developed  by  Preston  and  Koehler  (1968).  These  authors  considered  first  the 
mean-squared  angular  multiple-scattering  deflection  in  a path  length  s. 
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(14) 


<f)  = f%(s')ds' 

J O 

and  showed  that  the  corresponding  mean-squared  value  of  the  radial  deflection  is 

= f q(sO  (s-sO^  ds'  . 

J o 


The  function  q(s),  common  to  both  formulas,  was  evaluated  by  Preston  and  Koehler 
according  to  a prescription  given  by  Bethe  and  Ashkin  (1953).  Preston  and  Koehler  then 
approximated  the  radial  distribution  by  a Gaussian, 


f(p,s) 


P 


exp  (-p^  / <p^  ) , 


(16) 


and  used  this  Gaussian  as  a basis  for  calculating  central-axis  depth-dose  curves  for  circular  and 
for  rectangular  fields,  using  superposition  methods  equivalent  to  those  described  in  section  2. 

The  same  approach  was  also  used  by  Carlsson  and  Rosander  (1973),  who  calculated  radial 
absorbed-dose  distributions  in  water  from  a 185-MeV  proton  beam,  and  obtained  good  agreement 
with  radial  distributions  measured  by  them  with  a small  silicon  detector. 


Just  as  the  Gaussian  approximation  for  angular  multiple-scattering  deflections  can  be 
replaced  by  the  more  accurate  distribution  of  Moli^re  (1948),  the  Gaussian  ^proximation  for  the 
radial  distribution  f(p,s)  can  be  replaced  by  the  radial  distribution  of  Moli^re  (1955).  The 
evaluation  of  f(p,s)  thereby  becomes  somewhat  more  complicated,  but  can  still  be  done  very 
quickly  with  a computer.  The  relevant  equations  from  Moli^re’s  theory  are  given  in  seaion  4.2. 

Both  with  the  Gaussian  approximation  and  in  Moli^re’s  theory  the  energy  loss  of  the 
protons  can  be  taken  account  in  the  continuous-slowing-down  ^proximation:  at  every  point  along 
the  track  the  energy  loss  is  assumed  to  be  equal  to  the  stopping  power.  Energy-loss  straggling  is 
thus  neglected,  but  this  has  almost  no  effect  on  the  shape  of  the  radial  distribution.  This  was 
verified  by  comparative  calculations  with  PTRAN,  carried  out  with  and  without  energy-loss 
straggling.  As  illustrated  in  figure  5 for  the  case  of  a 160-MeV  beam,  energy  loss  straggling 
changes  the  radial  distribution  only  slightly  at  depths  z = 0.1  r^,  and  0.99  r^,,  but  practically  not 
at  all  at  depths  from  z = 0.2  r^  to  0.98  r^,. 

Moli^re’s  theory  provides  the  radial  distribution  f(p,s)  as  a function  of  the  path  length 
traversed  by  the  proton,  whereas,  strictly  speaking,  one  needs  f(p,z)  as  a function  of  the  depth  z. 
However,  the  difference  between  the  path  length  and  depth  is  exceedingly  small.  This  is 
illustrated  in  figure  6,  calculated  for  a 160-MeV  beam  with  a modified  version  of  PTRAN,  which 
shows  the  percentage  amount  by  which  the  average  path  length  s^v  differs  from  z. 

Direct  comparisons  between  radial  distributions  from  Moli^re’s  theory  and  histograms 
from  PTRAN  are  shown  in  figure  7a  (160  MeV)  and  in  figure  7b  (70  MeV),  at  depths 
z = 0.1  r^,  0.5  r^,  0.9  r^  and  0.99  r^.  The  agreement  is  generally  quite  close,  except  the  for  the 
depth  z = 0. 1 r^  where  a slight  difference  can  be  observed  near  the  peak  of  the  radial 
distribution.  The  origin  of  this  discrepancy  is  not  well  understood.  It  is  suspected  to  be  caused 
by  an  ^proximation  made  in  the  PTRAN  program  in  the  evaluation  of  later^  multiple-scattering 
displacements  [eqs  (2.7)  and  (2.8)  in  Berger,  1993a]. 
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One  the  basis  of  the  foregoing  considerations  and  comparison,  one  can  conclude  that 
Moli^re’s  theory  is  the  preferable  tool  for  calculating  the  radial  distribution  f(p,z),  and  should  be 
supplemented  by  the  more  costly  Monte  Carlo  method  only  to  obtain  results  at  depths  greater 
than  0.98  r^.  It  is  expected,  however,  that  the  Monte  Carlo  method  will  retain  its  usefulness 
even  at  shallow  and  intermediate  depths,  for  the  calculation  of  proton  transport  in  inhomogeneous 
media. 

4.2  Radial  Distribution  from  Molibre’s  Theory 

The  equations  are  listed  here  which  are  used  in  program  MORAD  to  calculate  the  radial 
distribution  f(p,s)  according  to  the  theory  of  Molibre.  MORAD  is  further  described  in 
section  6.3. 


Consider  a proton  starting  out  along  the  z axis.  Let  6 be  the  angular  multiple-scattering 
deflection,  i.e.,  the  polar  angle  specifying  the  direction  of  motion  after  the  traversal  of  a path 
length  s.  Molibre  (1948)  introduced  a scaled  angular  variable 


= 


d 

X^^/B  ’ 


(17) 


where  Xc  and  y/B  depend  on  the  path  length  s and  the  proton  energy.  In  terms  of  this  scaled 
angle,  the  Molibre  multiple-scattering  distribution  is 


= t?dtJ  + l + ...  i 

I B b2  J 

with  expansion  coefficients 

f ^ j “ ydy  exp(-y2/4) 


n log  r 

4^4 


(18) 


(19) 


where  denotes  a Bessel  function. 


Let  p be  the  radial  multiple-scattering  deflection,  i.e.,  the  distance  from  the  z-axis  which 
the  particle  has  reached  after  traversing  a path  length  s.  It  was  shown  by  Molibre  (1955),  that 
the  distribution  of  p is  given  by  the  scaled  distribution,  eq  (18),  provided  one  sets 


= 


p/s 


(20) 


Here  x^,  and  are  quantities  related  to,  but  slightly  different  from  the  corresponding  quantities 
Xc  and  B used  in  the  distribution  of  angular  multiple-scattering  deflections. 


*Actually  Moli^re  showed  that  the  distribution  of  x/s  (where  x is  the  lateral  displacement  in  the 
x-direction)  has  the  same  functional  form  as  the  distribution  of  the  scaled  projected  multiple  scattering 
deflection.  An  analogous  relation  holds  between  the  distribution  of  p/s  and  the  distribution  of  the  scaled 
spatial  multiple-scattering  deflection. 
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The  evaluation  of  and  must  be  done  for  path  lengths  s in  which  protons  lose  a 
considerable  fraaion  of  their  ener^.  In  Moli^re’s  theory  this  energy  loss  can  be  taken  account 
in  the  continuous-slowing-down  approximation.  The  proton  energy  is  expressed  as  a function  of 
the  path  length  traveled,  according  to  the  relation 


s = f 


ds' 


T - dE/dx 


(21) 


where  T^,  is  the  initial  energy,  Tg  is  the  energy  after  the  traversal  of  a path  length  s,  and  where 
-dE/dx  is  the  stopping  power. 


For  a compound 

2 2 
Xc  = E Wj  X,j  , 

j 

where  Wj  is  the  fraction  by  weight  of  the  atomic  constituent,  and  where 


= i:  h(s') 


l-s' 


ds'  , 


and 


h(s)  = 47r  N, 


m T+1 


M t(t+2)_ 


Aj 


(22) 


(23) 


(24) 


is  the  Avogadro  constant,  Zj  and  Aj  are  the  atomic  number  and  weight  of  the  constituent, 
r^  is  the  classical  electron  radius,  m/M  is  the  electron-proton  mass  ratio,  and  r is  the  proton 
kinetic  energy  in  units  of  the  proton  rest  mass  Mc^. 


The  parameter  is  obtained  as  the  solution  of  the  equation 
Bp  - log  Bp  = log(xp/xf)  + 1+27  , 


(25) 


where  7 = 0.5772156649...  is  Euler’s  constant,  and  the  screening  parameter  x is  given  by 


>°g  ’‘a  = 4 S jj  h(s')  " [log  Gj(s')  - Fj(s')/Zj]  ds'  , 


Xc  J 


1-s' 


(26) 


where 


M 0.88534 


.2/3 


1.13  + 3.76  (ZjO/^y 


(27) 


t(t  + 2) 


and  a is  the  fme-structure  constant  and  jS  is  the  proton  speed  in  units  of  the  speed  of  light. 

The  factor  kgp,  which  is  a function  of  Zj  a//5,  converts  Moli^re’s  result,  obtained  with  a 
Thomas-Fermi  potential,  to  a corresponding  result  for  a Hartree-Fock  potential.  Values  of  kjjp 
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are  given  in  Berger  (1993a).  Those  for  hydrogen  are  smaller  than  unity  by  up  to  25  percent,  and 
those  for  oxygen  are  larger  than  unity  by  up  to  12  percent. 

The  term  Fj/Zj  in  eq  (26)  is  a correction  due  to  Fano  (1956)  that  takes  into  account  the 
influence  of  the  electrons  of  the  target  atoms.  Fj  is  given  by 

Fj  = log  [ll30  (1  -Uj  - iP-n.  , (28) 

where  the  constant  Uj  has  the  value  —3.6  for  hydrogen  and  —5.1  for  oxygen. 

The  factor  (l-s7s)^  in  eqs  (23)  and  (26)  is  the  analog  of  the  factor  (s-s')^  in  eq  (15).  If 
this  factor  were  omitted,  one  would  obtain  instead  of  x^.  and  Bp  the  corresponding  quantities  Xc 
and  B for  the  distribution  of  angular  multiple-scattering  deflections. 


Table  2 lists  values  of  B^,  and  table  3 values  of  x^^^,  as  functions  of  z/r^,  for  seven 
energies  T^  between  250  MeV  and  50  MeV.  It  can  be  seen  that  these  quantities  are  very  slowly 
varying  functions  of  T^. 


The  width  of  the  Moli^re  distribution  can  be  characterized  as  a reduced  angle  1?^/^  at 
which  the  distribution  has  fallen  to  1/e  of  its  maximum.  Hanson  et  al.  (1951)  found  that  can 
be  approximated  as 


- (1.2/B^)  . 


(29) 


A more  accurate  expression  is 

0.80209  - 0.58365  B.  + 0.44997  hj 

_ P P 

^i/e  - 5 

1 - 0.32791  Bp  + 0.45000  B^ 


(30) 


Equations  (29)  and  (30)  can  of  course  also  be  applied  to  B^.  The  value  of  from  eq  (30)  is 
smaller  than  that  from  eq  (29)  by  1.49  percent  for  B^  = 4.5,  and  by  0.28  percent  for  B^  = 10, 
and  is  larger  by  0.43  percent  for  B^  = 20.  The  Gaussian  approximation  to  the  Molibre 
distribution  as  function  of  the  reduced  angle  is 


2x  fMG('’)  = exp(-t)2/4e) 


(31) 


In  figure  8,  the  Molibre  distribution  is  compared  with  the  Gaussian  approximation  for 
four  values  of  B^  (or  B)  between  20  and  4.5  (the  lowest  value  for  which  Moli^re’s  theory 
remains  accurate).  The  peak  value  of  the  Gaussian  approximation  is  too  high  and  its  tail  too  low. 
The  greater  the  value  of  B^  (or  B),  the  better  is  the  Gaussian  approximation. 

The  radial  deflection  py^  where  the  radial  distribution  has  fallen  to  1/e  of  its  maximum  is 

Pl/e  = S ^ <^l/e  • 
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Values  of  in  water  were  measured  by  Preston  and  Koehler  (1968).  Table  4 lists  these 
results,  together  with  the  theoretical  values  calculated  by  these  authors  in  the  Gaussian 
approximation,  and  with  values  from  the  Moli^re  theory.  A substantial  correction  was  applied  by 
Preston  and  Koehler  to  take  into  account  the  radial  spread  of  the  beam  incident  on  the  scattering 
target.  With  this  correction,  there  is  good  agreement  between  the  measured  and  calculated 
results. 

4.3  Database  of  Radial  Distributions  from  PTRAN 

The  results  of  calculations  with  PTRAN,  for  energies  T^  = 250,  200,  160,  130,  100,  70 
and  50  MeV,  obtained  with  a sample  of  1 million  proton  histories  in  each  case,  were  used  to 
prepare  a database  of  the  scaled  radial  distribution  2xp  zf(p,z)  as  function  of  p/z.  This  was  done 
for  depths  at  depths  at  which  Moli^re’s  theory  can  no  longer  be  used.  The  database  is  used  by  a 
program  called  PTRAD  which  generates,  by  linear  interpolation  with  respect  to  logT^,  radial 
distributions  f(p,z)  at  scaled  depths  z/r^  = 0.985,  0.99,  0.995,  1.0,  1.005,  1.01,  1.015,  1.020 
and  1.025  r^.  More  information  about  PTRAD  is  given  in  section  6.4. 


5.  Reduction  Factors  and  Absorbed-Dose  Distributions 

This  section  contains  illustrative  results  pertaining  to  an  incident  160-MeV  proton  beam 
(r^  = 17.65  cm  of  water).  The  reduaion  factors  and  absorbed-dose  distributions  were  obtained 
with  depth-dose  distributions  dD/dz  from  program  PTPOL,  and  with  radial  distributions  f(p,z) 
from  programs  MORAD  and  PTRAD.  The  reduction  factors  for  circular  fields  were  generated 
with  program  FCIR,  and  those  for  rectangular  fields  with  program  FREC.  Unless  the  contrary  is 
stated,  the  density  of  pencil  beams  is  assumed  to  constant  over  the  field  area.  Absorbed-dose 
values  are  given  in  units  of  MeV/g,  normalized  to  one  incident  proton  per  cm^. 

5.1  Circular  Fields 

Figures  9a-d  show  the  reduction  factors  [defined  by  eqs  (5a,b)]  at  depths  z = 0.1  r^, 

0.5  r^  and  0.99  r^,  as  funaions  of  the  distance  r from  the  center  of  circular  fields  with  radii 
R = 8,  4,  2 and  1 mm,  respectively.  Figure  10  shows  the  reduction  faaors  as  function  of  depth 
along  the  central  axis,  again  for  field  radii  R = 8,  4,  2,  and  1 mm. 

Central-axis  absorbed-dose  distributions  for  field  radii  R = 4,  2 and  1 mm  are  plotted  in 
figure  11.  Also  plotted  for  comparison  is  the  absorbed  dose  from  an  unbounded  field.  This 
figure  shows  how  the  reduction  of  the  field  radius  reduces  the  magnitude  of  the  Bragg  peak,  in 
agreement  with  earlier  findings  of  Preston  and  Koehler  (1968). 

Figures  12a,b,c  shows  absorbed  dose  distributions  as  functions  of  depth,  along  lines 
parallel  to  the  z axis,  for  various  distances  r between  these  lines  and  the  z axis,  for  circular  fields 
with  radii  R = 4,  2 and  1 mm,  respectively. 

Figure  13  shows  the  average  absorbed  dose  as  a function  of  depth,  for  circular  fields  with 
radius  R = 8,  4,  2 and  1 mm.  The  absorbed  dose  is  averaged  over  thin  disks  of  radius  R 
centered  on  the  z axis. 


11 


Figures  14a,b,c  show  reduction  factors  at  depths  z = 0.1  r^,  0.5  r^,,  and  0.9  r^  for  a 
circular  field  (radius  R = 4 mm)  calculated  for  a field  density  that  is  not  constant  but  is  given  by 
n(r)  = exp(-Qr^),  where  r is  the  distance  from  the  center  of  the  field.  These  results  were 
calculated  from  eqs  (7a,b),  for  parameter  values  Q = 2,  1,  and  0 (constant  density). 

5.2  Rectangular  Fields 

Figures  15a,b,c  compare  reduction  factors  for  circular  fields  with  reduction  factors  for 
square  fields,  at  depths  z = 0.1  r^,  0.5  r^  and  0.9  r^,  respectively.  The  circular  fields  and 
square  fields  have  die  same  area  (64  mm^).  Two  sets  of  curves  are  shown  for  the  square  field. 

In  one  case  (indicated  as  0 deg),  the  reduction  factors  are  plotted  along  a line  that  starts  at  the 
center  of  the  field  and  is  parallel  to  an  edge  of  the  square.  In  the  other  case  (indicated  as  45  deg) 
the  line  that  starts  at  the  center  of  the  field  and  passes  though  a comer  of  the  square. 

Figure  16  shows  the  absorbed  dose  along  the  central  axis,  for  square  fields  with  sides 
equal  to  8,  4,  2 or  1 mm.  Figure  17  shows  similar  results  for  a field  in  the  form  of  a slit,  i.e., 
with  one  side  of  the  rectangular  field  effectively  unbounded,  and  with  the  other  side  having  a 
width  equal  to  8,  4,  2 or  1 mm. 

Figures  18a,b,c  show  reduction  factors  and  absorbed-dose  distributions  for  a rectangular 
1 mm  X 10  mm  field.  Figure  18a  shows  this  quantities  as  a function  of  depth,  along  the  central 
field  axis.  Figure  18b  shown  them  at  a depth  z = 0.5  r^,,  as  a function  of  radial  distance  r from 
the  field  axis,  where  the  radial  distance  is  measured  along  a line  that  goes  through  the  center,  and 
is  parallel  to  the  long  side,  of  the  rectangle.  Figure  18c  is  similar  to  18b,  except  that  the  radial 
distance  is  measured  along  a line  that  is  parallel  to  the  sort  side  of  the  rectangle. 


6.  Computer  Programs 

The  program  files,  data  file  and  sample  output  files  supplied  with  this  report  are  stored  on 
a single  3.5"  1.44-Mb  disk.  The  program  files  include  Fortran  source  code  as  well  as  executable 
code  for  use  with  an  IBM-compatible  personal  computer.  The  source  code,  data  files,  and 
sample  output  files  are  written  in  ASCII  format,  and  can  be  transferred  to  other  types  of 
computers.  A listing  of  all  files  can  be  found  in  Appendbc  1,  and  the  source  code  for  all  the 
Fortran  programs  in  Appendix  2. 

The  programs  prompt  the  user  to  supply  input  data  and/or  to  indicate  the  names  of  files 
which  contain  the  required  input  data.  The  output  files  are  supplied  with  headings  and 
descriptions,  so  that  they  are  practically  self-explanatory. 

6.1  Program  PTPOL 

PTPOL  calculates  the  energy  deposition  distribution  dD/dz  in  water  for  a monoenergetic 
proton  beam.  The  user  must  supply  the  beam  energy  (in  MeV)  and  the  name  of  the  output  file. 

If  the  desired  energy  is  smaller  than  50  MeV  or  greater  than  250  MeV,  the  program  halts  and 
tells  the  user  that  the  energy  is  out  of  bounds.  The  output  file  includes  a list  of  depths,  the 
values  of  dD/dz  at  these  depths,  in  units  of  MeV/g,  and  the  corresponding  relative  depth-dose 
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values  (with  unit  peak  height).  Also  included  in  the  ou^ut  are  the  depth  at  which  the  peak 
occurs,  and  the  peak  value  of  dD/dz.  The  sample  ouq)ut  file  PTPOL.160,  shown  in  ti)le  5,  is 
for  a beam  energy  of  160  MeV. 

6.2  Program  PTGPOL 

PTGPOL  calculates  the  energy  deposition  distribution  dD/dz  in  water  for  a proton  beam 
with  a Gaussian  energy  spectrum.  The  user  must  specify  the  average  energy  T^^  and  the 
standard  deviation  Oj  of  die  spectrum.  The  latter  must  be  entered  in  terms  of  a relative  standard 
deviation  P = 100  Tlie  Gaussian  distribution  is  truncated  at  the  energies 

Tj  2 = T^v  ± ^ ^T-  cut-off  parameter  q (recommended  value  4)  must  be  specified  by  the 
user.  If  either  T^  or  T2  are  smaller  than  50  MeV  or  greater  than  250  MeV,  the  program  halts 
and  informs  the  user  that  these  energies  are  out  of  bounds. 

The  set  of  depths  at  which  dD/dz  is  to  be  evaluated  must  be  supplied  in  a previously- 
prepared  file.  The  first  entry  in  this  file  must  be  the  number  of  depths  (no  greater  than  201),  and 
the  other  entries  must  be  depth-values  expressed  in  units  of  the  CSDA  range  at  energy  T^^.  A 
default  file  ZRLIST,  supplied  with  the  program,  has  52  depth-values. 

The  output  file  from  PTGPOL  contains  similar  information  as  the  output  file  from 
PTPOL,  and  in  addition  lists  the  range  at  energy  T^^,  and  the  parameters  of  the  Gaussian  beam 
spectrum.  A sample  output  file  PTGPOL.OUT,  shown  in  table  6,  pertains  to  a beam  with  a 
mean  energy  of  160  MeV  and  a relative  standard  deviation  of  P = 1 percent. 

6.3  Program  MORAD 

MORAD  calculates  the  radial  distribution  f(p,s)  according  to  the  theory  of  Moli^re,  using 
the  equations  listed  in  section  4.2.  In  response  to  prompts,  the  user  must  specify  a)  the  beam 
energy  (in  MeV),  b)  a set  of  path  lengths,  and  c)  the  set  of  radial  distances  at  which  the  radial 
distributions  are  to  be  evaluated  for  each  path  length. 

The  path  lengths  (in  units  of  the  CSDA  range)  can  be  read  from  a previously  prepared 
file.  The  first  line  of  such  file  should  contain  NMAX,  the  number  of  path  lengths,  and 
subsequent  lines  should  contain  the  path  lengths,  separated  by  blanks  or  other  delimeter. 
Alternatively,  the  path  lengths  can  be  specified  in  terms  of  a maximum  value  ZRMAX  and  a 
number  NMAX,  from  which  MORAD  computes  the  path  lengths  n(ZRMAX/NMAX),  n = 1,  2, 
...,  NMAX.  It  is  also  possible  to  specify  that  ZRMAX  should  have  the  default  value  0.98,  and 
NMAX  the  default  value  98. 

The  Moli^re  distribution,  as  a fimction  of  the  reduced  angle  t?,  extends  from  zero  to 
infinity.  This  is  unrealistic  and  is  due  to  an  ^proximation  made  in  the  theory.  In  MORAD  the 
distribution  is  truncated  at  renormalized  to  unity.  Prior  to  renormaliza- 

tion, the  integral  over  the  distribution  up  to  is  smaller  than  unity  by  a fraction  of  a percent, 
and  is  printed  out  as  part  of  the  ouq)ut  of  MORAD.  The  corresponding  truncation  value  of  the 
radial  defiection  is 
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(33) 


^max  ” S 


max 


The  radial  distribution  is  tabulated  by  MORAD  at  KMAX  equidistant  value  of  between  0 
and  Pmax-  KMAX  is  an  input  parameter.  Table  7 shows  an  excerpt  from  a sample  output  file 
MORAD.  160  for  98  depths  (path  lengths),  and  table  8 an  ou4)ut  file  MORAD160.050  for  a 
single  depth  (z  = 0.5  r^). 

6.4  Program  PTRAD 

PTRAD  calculates  the  radial  distribution  f(p,z)  using  Monte  Carlo  results  from  PTRAN 
stored  in  a database  prepared  as  described  in  section  4.3.  There  are  9 files  in  this  database, 
called  RADF.n,  n = 32,  33,. ..,40,  which  contain  scaled  radial  distributions  at  9 depths 
(z  = 0.985,  0.99,  0.995,  1.0,  1.005,  1.01,  1.015,  1.020  and  1.025  r^).  The  distributions  are 
generated  at  201  equidistant  values  of  p which  are  listed  in  the  output  file.  In  order  to  run 
PTRAD,  the  user  merely  has  to  specify  the  beam  energy,  in  MeV,  and  the  name  of  the  output 
file.  A sample  output  file  PTRAD.  160  is  shown  in  table  9. 

6.5  Program  FCIR 

FCIR  uses  the  equations  given  in  section  2.1  to  calculate  reduction  factors  and  absorbed- 
dose  values  as  functions  of  the  depth  z and  of  the  radial  distance  r from  the  center  of  a circular 
field.  The  user  must  specify  the  field  radius,  in  cm,  and  the  set  of  distances  r for  which  tiie 
calculation  is  to  be  done.  The  radial  distances  can  be  read  from  a prepared  file,  entered  from  the 
keyboard,  or  specified  in  terms  of  a maximum  radial  distance  r^^^  and  a number  LMAX  of 
radial  distances.  With  the  third  choice,  FCIR  computes  a set  of  I^AX+1  distances 

f = 0,  1,  2,  ...,  LMAX.  If  LMAX  is  1,  calculations  are  made  only  for 
r = 0,  along  the  central  axis  of  the  field. 

The  user  is  queried  whether  the  requested  path  lengths  are  all  smaller  than  0.98  r^,.  If  the 
answer  is  yes,  input  has  to  be  supplied  only  from  a file  generated  by  MORAD.  If  the  answer  is 
no,  input  must  in  addition  be  supplied  from  a file  generated  by  PTRAD.  The  beam  energy  and 
path  lengths,  and  the  radial  distributions  for  pencil  beams,  are  supplied  to  FCIR  via  the  files  from 
MORAD  and  PTRAD.  Finally,  the  user  must  supply  the  name  of  the  output  file. 

Sample  output  files  FCIR.l  and  FCIR.2  , for  a circular  field  with  a radius  of  2 mm,  are 
shown  in  tables  10  and  11,  respectively.  FCIR.l  contains  reduction  factors  and  absorbed-dose 
values,  as  a function  of  depth,  along  the  central  axis.  FCIR.2  contains  reduction  factors  and 
absorbed-dose  values  as  a function  of  the  radial  distance  r,  at  a depth  z = 0.5  r^. 

6.6  Programs  FREC,  FRECl,  and  FREC2 

FREC  uses  the  equations  given  in  section  2.2  to  calculate  reduction  factors  and  absorbed- 
dose  values  as  functions  of  the  depth  z and  of  the  distance  r from  the  center  of  a rectangular 
field.  The  field  is  assumed  to  extend  from  -A  to  A in  x,  and  from  -B  to  B in  y,  with  the  origin 
at  the  center  of  the  rectangular  field.  The  user  must  specify  the  parameters  A and  B,  and  also 
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the  angle  ot  with  respect  to  the  x-axis  of  a line  starting  at  the  center  of  the  field  along  wiiich 
reduaion  factors  and  absorbed-dose  values  are  to  be  calculated.  The  user  can  also  indicate  that 
the  line  should  pass  through  a comer  of  the  field,  in  vMch  case  FREC  computes  a. 

The  remaining  input  data,  including  the  radial  distances  r along  the  chosen  line,  and  the 
input  files  from  MORAD  and  PTRAD,  with  beam  energy,  path  lengths  and  radial  distributions 
from  pencil  beams,  must  be  supplied  in  the  same  manner  as  described  in  section  6.5  for  FCIR. 

Sample  ou4)ut  files  FREC.l,  FREC.2A  and  FREC.2B  are  shown  in  tables  12,  13,  and 
14,  which  pertain  to  a rectangular  1 mm  x 10  mm  field.  FREC.l  contains  reduction  factors  and 
absorbed-dose  values  as  a function  of  depth,  along  the  central  axis.  FREC.2A  and  FREC.2B 
contain  reduction  factors  and  absorbed-dose  values  at  a depth  0.5  r^,  along  lines  that  pass  through 
the  center  of  the  field  and  are  parallel  to  the  short,  respectively,  the  long  side  of  the  rectangular 
field. 


Program  FRECl  is  similar  to  FREC,  but  provides  reduction  factors  and  absorbed-dose 
values  not  along  a specified  line,  but  for  a set  of  specified  points  (x,y).  Table  15  shows  a 
sample  output  file  for  a rectangular  5 mm  x 2 mm  field. 

Program  FREC2  performs  the  same  function  as  FRECl,  but  evaluates  eq  (1)  with  a 
double  numerical  quadrature,  with  respect  to  x from  -A  to  A,  and  with  respect  to  y from  y from 
-B  to  B.  The  integration  limits  are  specified  by  the  external  functions  FUNl  and  FUN2.  The 
pencil-beam  density,  n(x,y)  = 1,  is  specified  by  the  external  function  FUN3.  Program  FREC2 
could  easily  be  modified  so  as  to  apply  to  a field  of  different  size  and  sh^e,  by  changing  the 
functions  FUNl  and  FUN2.  Any  pencil-beam  density  n(x,y)  could  be  used  by  the  choice  of  an 
appropriate  function  FUN3. 

6.7  Programs  AXPLOT  and  RADPLOT 

AXPLOT  plots  reduaion  faaors  and  absorbed-dose  values  along  the  central  axis  of  a 
circular  or  rectangular  field,  using  input  files  generated  by  FCIR  or  F^C.  RADPLOT  plots 
reduaion  factors  and  absorbed-dose  values  at  a fixed  depth,  as  functions  of  the  radial  distance 
from  the  center  of  a circular  or  rectangular  field.  These  programs  must  be  compiled  and  linked 
with  a library  from  a graphics  program  by  Kahaner  and  Anderson  (1990).  Figure  18  of  this 
report  was  produced  with  AXPLOT,  and  figures  19  and  20  with  RADPLOT. 
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Table  1.  Stopping  powers  and  CSDA  ranges  for  protons  in  water.  Calculated  with  program 
PSTAR  (Berger,  1992). 


T dE/dx  r^ 

(MeV)  (MeV  cm^/g)  g/cm^ 


250.0 

237.5 

225.0 

212.5 

200.0 

190.0 

180.0 

170.0 

160.0 

152.5 

145.0 

137.5 

130.0 

122.5 

115.0 

107.5 

100.0 

92.5 

85.0 

77.5 

70.0 

65.0 

60.0 

55.0 

50.0 


3.910 

4.032 

4.169 

4.320 

4.491 

4.642 

4.810 

4.997 

5.207 

5.381 

5.573 
5.784 
6.018 
6.280 

6.574 

6.907 

7.286 

7.723 

8.233 

8.834 

9.555 

10.121 

10.775 

11.539 

12.443 


37.94 

34.79 
31.74 

28.80 
25.96 

23.77 

21.65 

19.614 

17.653 

16.237 

14.867 

13.546 

12.275 

11.055 

9.888 

8.775 

7.718 

6.718 
5.777 
4.897 

4.080 
3.572 
3.093 
2.644 
2.227 
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Table  2.  Parameter  in  Moli^re’s  radial  distribution  (for  water),  as  a function  of  s/r^,  where  s is 
the  path  len^  and  r^  is  the  CSDA  range  path  length,  for  various  values  of  die  initial 
proton  energy  T^,. 


(MeV)  250.0 

200.0 

160.0 

130.0 

100.0 

70.0 

50.0 

r^  (g/cm^)  37.94 

25.95 

17.65 

12.28 

7.718 

4.080 

2.227 

s/fo 


0.00006 

6.284 

5.979 

5.674 

5.391 

5.032 

4.536 

4.055 

0.00008 

6.624 

6.322 

6.022 

5.742 

5.388 

4.901 

4.431 

0.00010 

6.886 

6.586 

6.288 

6.011 

5.660 

5.180 

4.717 

0.00015 

7.358 

7.061 

6.767 

6.494 

6.148 

5.677 

5.225 

0.00020 

7.690 

7.395 

7.103 

6.832 

6.490 

6.024 

5.578 

0.00030 

8.154 

7.862 

7.573 

7.304 

6.967 

6.506 

6.067 

0.0004 

8.481 

8.191 

7.903 

7.637 

7.301 

6.845 

6.410 

0.0005 

8.733 

8.444 

8.158 

7.893 

7.559 

7.105 

6.674 

0.0006 

8.939 

8.651 

8.365 

8.101 

7.769 

7.317 

6.887 

0.0008 

9.262 

8.975 

8.691 

8.429 

8.098 

7.649 

7.223 

0.0010 

9.512 

9.226 

8.943 

8.681 

8.352 

7.905 

7.481 

0.0015 

9.964 

9.680 

9.398 

9.138 

8.811 

8.368 

7.947 

0.002 

10.283 

10.000 

9.720 

9.461 

9.135 

8.694 

8.275 

0.003 

10.732 

10.450 

10.171 

9.913 

9.589 

9.150 

8.735 

0.004 

11.049 

10.768 

10.489 

10.233 

9.910 

9.473 

9.059 

0.005 

11.294 

11.014 

10.736 

10.480 

10.158 

9.722 

9.310 

0.006 

11.494 

11.214 

10.937 

10.681 

10.360 

9.926 

9.514 

0.008 

11.809 

11.530 

11.254 

10.999 

10.679 

10.245 

9.836 

0.010 

12.053 

11.775 

11.499 

11.245 

10.925 

10.493 

10.084 

0.015 

12.496 

12.218 

11.943 

11.690 

11.372 

10.941 

10.535 

0.020 

12.809 

12.532 

12.258 

12.006 

11.688 

11.259 

10.853 

0.030 

13.251 

12.975 

12.702 

12.450 

12.133 

11.706 

11.302 

0.040 

13.564 

13.289 

13.016 

12.765 

-'12.449 

12.022 

11.619 

0.050 

13.807 

13.532 

13.260 

13.009 

12.694 

12.268 

11.866 

0.06 

14.006 

13.731 

13.460 

13.209 

12.895 

12.469 

12.068 

0.08 

14.320 

14.046 

13.775 

13.525 

13.211 

12.787 

12.386 

0.10 

14.565 

14.291 

14.021 

13.771 

13.458 

13.034 

12.634 

0.15 

15.013 

14.740 

14.470 

14.221 

13.909 

13.486 

13.088 

0.20 

15.334 

15.061 

14.792 

14.544 

14.232 

13.811 

13.413 

0.30 

15.795 

15.524 

15.256 

15.008 

14.698 

14.278 

13.881 

0.4 

16.132 

15.862 

15.595 

15.348 

15.039 

14.620 

14.224 

0.5 

16.404 

16.135 

15.868 

15.623 

15.314 

14.896 

14.501 

0.6 

16.636 

16.368 

16.103 

15.858 

15.550 

15.133 

14.738 

0.8 

17.040 

16.774 

16.511 

16.267 

15.961 

15.546 

15.152 

0.9 

17.234 

16.970 

16.708 

16.466 

16.160 

15.746 

15.352 

0.99999 

17.461 

17.199 

16.939 

16.698 

16.394 

15.980 

15.587 
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Table  3.  Scale  parameter  x,.  in  Moli^re’s  radial  distribution  (for  water)  as  function  of  s/r^, 

where  s is  the  path  length  and  r^  is  the  CSDA  range,  for  various  values  of  the  initial  proton 
energy  T^. 


T„(MeV) 

250.0 

200.0 

160.0 

130.0 

100.0 

70.0 

50.0 

To  (g/cm^) 

37.94 

25.95 

17.65 

12.28 

7.718 

4.080 

2.227 

S/Fo 

0.00006 

1.511 

1.495 

1.477 

1.458 

1.433 

1.394 

1.350 

0.00008 

1.551 

1.537 

1.521 

1.505 

1.483 

1.449 

1.411 

0.00010 

1.581 

1.569 

1.554 

1.540 

1.520 

1.489 

1.456 

0.00015 

1.635 

1.624 

1.613 

1.601 

1.584 

1.559 

1.532 

0.00020 

1.671 

1.662 

1.652 

1.642 

1.627 

1.606 

1.583 

0.00030 

1.721 

1.714 

1.706 

1.698 

1.686 

1.669 

1.651 

0.0004 

1.755 

1.750 

1.743 

1.736 

1.726 

1.712 

1.697 

0.0005 

1.781 

1.776 

1.771 

1.765 

1.756 

1.744 

1.732 

0.0006 

1.802 

1.798 

1.793 

1.788 

1.781 

1.770 

1.759 

0.0008 

1.834 

1.832 

1.828 

1.824 

1.818 

1.810 

1.801 

0.0010 

1.859 

1.857 

1.854 

1.851 

1.846 

1.840 

1.833 

0.0015 

1.903 

1.902 

1.901 

1.899 

1.897 

1.893 

1.890 

0.0020 

1.933 

1.934 

1.933 

1.932 

1.931 

1.930 

1.929 

0.0030 

1.975 

1.977 

1.978 

1.978 

1.979 

1.980 

1.982 

0.0040 

2.004 

2.007 

2.009 

2.010 

2.012 

2.015 

2.018 

0.0050 

2.027 

2.030 

2.033 

2.035 

2.037 

2.042 

2.046 

0.0060 

2.045 

2.049 

2.052 

2.054 

2.058 

2.063 

2.069 

0.0080 

2.073 

2.078 

2.082 

2.085 

2.090 

2.097 

2.104 

0.010 

2.095 

2.100 

2.105 

2.109 

2.114 

2.122 

2.131 

0.010 

2.135 

2.141 

2.147 

2.152 

2.159 

2.169 

2.180 

0.020 

2.163 

2.170 

2.176 

2.182 

2.190 

2.201 

2.214 

0.030 

2.203 

2.211 

2.218 

2.225 

2.234 

2.248 

2.262 

0.040 

2.232 

2.240 

2.249 

2.256 

2.266 

2.281 

2.297 

0.050 

2.255 

2.264 

2.273 

2.281 

2.292 

2.307 

2.325 

0.06 

2.274 

2.284 

2.293 

2.302 

2.313 

2.330 

2.348 

0.08 

2.306 

2.316 

2.326 

2.335 

2.348 

2.366 

2.385 

0.10 

2.332 

2.343 

2.354 

2.363 

2.376 

2.395 

2.416 

0.15 

2.385 

2.397 

2.408 

2.419 

2.433 

2.454 

2.477 

0.20 

2.428 

2.441 

2.453 

2.464 

2.480 

2.502 

2.526 

0.30 

2.503 

2.517 

2.530 

2.543 

2.560 

2.585 

2.611 

0.4 

2.573 

2.588 

2.602 

2.615 

2.633 

2.660 

2.688 

0.5 

2.643 

2.658 

2.673 

2.688 

2.707 

2.735 

2.765 

0.6 

2.716 

2.732 

2.748 

2.763 

2.783 

2.813 

2.845 

0.8 

2.885 

2.903 

2.921 

2.938 

2.960 

2.994 

3.030 

0.9 

2.993 

3.012 

3.031 

3.049 

3.073 

3.110 

3.148 

0.99999 

3.143 

3.163 

3.184 

3.204 

3.230 

3.271 

3.313 
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Table  4.  Comparison  of  the  calculated  radial  spread  of  proton  beams  in  water  with  experimental 
results  of  Preston  and  Koehler  (1968).  The  radial  spread  is  defined  as  the  radial  distance 
from  the  beam  axis  at  which  intensity  has  fallen  to  a fraction  1/e  of  the  peak  intensity. 


Beam  Energy 

(MeV) 

Depth 

(cm) 

Pm 

Pa 

Radial  Spread  (mm), 

Pc 

PPK 

Pm 

127 

5.7 

1.90 

1.51 

1.15  ± 0.15 

1.08 

1.02 

127 

8.7 

2.70 

1.60 

2.18  ± 0.10 

2.17 

2.07 

127 

11.4 

3.85 

1.68 

3.46  ± 0.09 

3.50 

3.40 

134 

12.4 

3.80 

0.94 

3.68  ± 0.10 

3.69 

3.65 

Pc 


Pm 

Pa 

(Pm"  - Pa")'"" 
PPK 
Pm 


measured 

measured  in  absence  of  scattering  material 
corrected  experimental  value 

calculated  by  Preston  and  Koehler,  using  Gaussian  approximation 
calculated  from  Moli^re’s  theory. 
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Table  5.  Output  file  from  program  PTPOL. 


Program  PTPOL,  output  file  PTPOL. 160 

EBEAM  = average  beam  energy,  MeV 

RANGE  = CSDA  range  at  energy  EBEAM,  g/cm^ 

ZM  = depth  (in  units  of  RANGE)  at  which  energy-deposition 
DMAX  = energy  deposition  per  unit  depth  at  peak,  MeV  cmVg 
NMAX  = number  of  depths 

EBEAM  RANGE  ZM  DMAX  NMAX 

152.00  16.1474  0.9903  31.2816  44 

curve  peaks 

Depths,  in 

units  of  RANGE 

0.000 

0.100  0.200 

0.300 

0.400 

0.450 

0.500 

0.550  0.600 

0.650 

0.700 

0.720 

0.740 

0.760  0.780 

0.800 

0.820 

0.840 

0.860 

0.880  0.900 

0.910 

0.920 

0.930 

0.940 

0.950  0.955 

0.960 

0.965 

0.970 

0.975 

0.980  0.985 

0.990 

0.995 

1.000 

1.005 

1.010  1.015 

1.020 

1.025 

1.030 

1.035 

1.040 

dD/dz , energy- 

deposition  distribution. 

MeV  cmVgj 

per  incident  proton 

6.379316 

6.443711 

6.545821 

6.693469 

6.904040 

7.041533 

7.207235 

7.408860 

7.656943 

7.967148 

8.364057 

8.555535 

8.771009 

9.015185 

9.294683 

9.617826 

9.996713 

10.449075 

10.998257 

11.686963 

12.582426 

13.143499 

13.813137 

14.630116 

15.659860 

17.021864 

17.898077 

18.979285 

20.356806 

22.178840 

24.537248 

27.325106 

29.988105 

31.277747 

29.860523 

25.003223 

17.941608 

10.807292 

5.297353 

2.073649 

0.639564 

0.164671 

0.031402 

0.004889 

Relative  energy-deposition 

distribution 

(peak  value 

unity) 

0.203932 

0.205991 

0.209255 

0.213975 

0.220706 

0.225102 

0.230399 

0.236844 

0.244775 

0.254692 

0.267380 

0.273501 

0.280389 

0.288195 

0.297130 

0.307460 

0.319572 

0.334033 

0.351589 

0.373606 

0.402231 

0.420168 

0.441574 

0.467691 

0.500610 

0.544150 

0.572161 

0.606724 

0.650761 

0.709007 

0.784400 

0.873521 

0.958651 

0.999878 

0.954573 

0.799296 

0.573552 

0.001004 

0.345484 

0.000156 

0.169344 

0.066290 

0.020445 

0.005264 
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Table  6.  Output  file  from  program  PTGPOL. 


Program  PTGPOL,  output  file  PTGPOL. 160 

List  of  depths  from  file  ZRFILE 

EAV  = average  beam  energy,  MeV 

SIGMA  = standard  deviation  of  beam  spectrum,  MeV 

P = 100*SIGMA/EAV,  relative  standard  deviation 

CUT  ==  cut-off  for  Gaussian  beam  (in  units  of  SIGMA) 

RGREF  * CSDA  range  at  energy  EAV 

ZM  = depth  (in  units  of  RANGE)  at  which  energy— deposition  curve  peaks 
UMAX  = energy  deposition  per  unit  depth  at  peak,  MeV  cmVg 

NRMAX  = number  of  depths 

EAV  SIGMA  P CUT  RGREF  ZM  UMAX  NRMAX 

160.00  1.60  1.00  4.00  17.6500  0.9813  22.5506  52 

Depths,  in  units  of  RGREF 

0.000000 

0.100000 

0.200000 

0.300000 

0.400000 

0.450000 

0.500000 

0.550000 

0.600000 

0.650000 

0.700000 

0.720000 

0.740000 

0.760000 

0.780000 

0.800000 

0.820000 

0.840000 

0.860000 

0.880000 

0.900000 

0.910000 

0.920000 

0.930000 

0.940000 

0.950000 

0.955000 

0.960000 

0.965000 

0.970000 

0.975000 

0.980000 

0.985000 

0.990000 

0.995000 

1.000000 

1.005000 

1.010000 

1.015000 

1.020000 

1.025000 

1.030000 

1.035000 

1.040000 

1.045000 

1.050000 

1.055000 

1.060000 

1.065000 

1.070000 

1.075000 

1.080000 

dD/dz,  energy— deposition  distribution. 

MeV  cmVg, 

per  incident 

proton 

6.245952 

6.290157 

6.370409 

6.493944 

6.677337 

6.799634 

6.948656 

7.131837 

7.359206 

7.645930 

8.015577 

8.194778 

8.396866 

8.626770 

8.890688 

9.197047 

9.557862 

9.990335 

10.521089 

11.192994 

12.083185 

12.654453 

13.354994 

14.249679 

15.445681 

17.075772 

18.072508 

19.162411 

20.280644 

21.318145 

22.125494 

22.531537 

22.376124 

21.549951 

20.029351 

17.893136 

15.313623 

12.523026 

9.765560 

7.250373 

5.118912 

3.433616 

2.186644 

1.321339 

0.757280 

0.411389 

0.211664 

0.102985 

0.047249 

0.020323 

0.008121 

0.002974 

Relative  energy- 

-deposition 

distribution 

(peak  value 

unity) 

0.276974 

0.278935 

0.282493 

0.287972 

0.296104 

0.301527 

0.308136 

0.316259 

0.326341 

0.339056 

0.355448 

0.363394 

0.372356 

0.382551 

0.394254 

0.407840 

0.423840 

0.443018 

0.466554 

0.496349 

0.535824 

0.561157 

0.592222 

0.631897 

0.684933 

0.757219 

0.801419 

0.849750 

0.899338 

0.945345 

0.981147 

0.999153 

0.992261 

0.955625 

0.888194 

0.793465 

0.679077 

0.555329 

0.433050 

0.321515 

0.226996 

0.152262 

0.096966 

0.058594 

0.033581 

0.018243 

0.009386 

0.004567 

0.002095 

0.000901 

0.000360 

0.000132 
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Table  7.  Excerpt  from  output  file  MORAD.160  from  program  MORAD 


Program  MORAD,  output  file  MORAD. 160 
98  depths 

201  radial  distances 


TIN 

energy  of  incident  proton  beam.  MeV 

RGIN 

= 

CSDA  range  at  energy  TIN,  g/cm^ 

ZR 

= 

depth,  in  units  of  RGIN 

RMAX 

= 

largest  radial  distance 

SUM 

= 

cumulative  integral  of  Moliere  distribution  up  to 

TIN  RGIN  ZR  RMAX  SUM 

160.000000  17.653381  0.010000  0.003716  0.999108 


Radial  distances,  p,  cm 


O.OOOOOE+OO 

1.85799E-05 

3.71597E-05 

5.57396E-05 

7.43194E-05 

9.28993E-05 

1.11479E-04 

1.30059E-04 

1.A8639E-04 

1.67219E-04 

1.85799E-04 

2.04378E-04 

2.22958E-04 

2.41538E-04 

2.60118E-04 

2,78698E-04 

2.97278E-04 

3.15858E-04 

3.34437E-04 

3.53017E-04 

3.71597E-04 

3.90177E-04 

4.08757E-04 

4.27337E-04 

4.45917E-04 

4,64496E-04 

4.83076E-04 

5.01656E-04 

5.20236E-04 

5.38816E-04 

5.57396E-04 

5.75976E-04 

5.94555E-04 

6.13135E-04 

6.31715E-04 

6.50295E-04 

6.68875E-04 

6.87455E-04 

7.06035E-04 

7.24615E-04 

7.43194E-04 

7.61774E-04 

7.80354E-04 

7.98934E-04 

8.17514E-04 

8.36094E-04 

8.54674E-04 

8.73253E-04 

8.91833E-04 

9.10413E-04 

9.28993E-04 

9,47573E-04 

9.66153E-04 

9.84733E-04 

1.00331E-03 

1,02189E-03 

1.04047E-03 

l,05905E-03 

1.07763E-03 

1.09621E-03 

l,11479E-03 

1.13337E-03 

1.15195E-03 

1.17053E-03 

1.18911E-03 

1.20769E-03 

1.22627E-03 

1.24485E-03 

1.26343E-03 

1.28201E-03 

1.30059E-03 

1.31917E-03 

1.33775E-03 

1.35633E-03 

1.37491E-03 

1.39349E-03 

1.41207E-03 

1.43065E-03 

1.44923E-03 

1.46781E-03 

1.48639E-03 

1.50497E-03 

1.52355E-03 

1.54213E-03 

1.56071E-03 

1.57929E-03 

1.59787E-03 

1.61645E-03 

1.63503E-03 

1.65361E-03 

1.67219E-03 

1.69077E-03 

1.70935E-03 

1.72793E-03 

1.74651E-03 

1.76509E-03 

1.78367E-03 

1.80225E-03 

1.82083E-03 

1.83941E-03 

1.85799E-03 

1.87657E-03 

1.89515E-03 

1.91373E-03 

1.93231E-03 

1.95089E-03 

1.96947E-03 

1.98804E-03 

2,00662E-03 

2.02520E-03 

2.04378E-03 

2.06236E-03 

2.08094E-03 

2.09952E-03 

2,11810E-03 

2.13668E-03 

2.15526E-03 

2.17384E-03 

2.19242E-03 

2.21100E-03 

2.22958E-03 

2.24816E-03 

2.26674E-03 

2.28532E-03 

2.30390E-03 

2.32248E-03 

2.34106E-03 

2.35964E-03 

2.37822E-03 

2.39680E-03 

2.41538E-03 

2.43396E-03 

2.45254E-03 

2.47112E-03 

2.48970E-03 

2.50828E-03 

2.52686E-03 

2.54544E-03 

2.56402E-03 

2.58260E-03 

2.60118E-03 

2.61976E-03 

2.63834E-03 

2.65692E-03 

2.67550E-03 

2.69408E-03 

2.71266E-03 

2.73124E-03 

2.74982E-03 

2.76840E-03 

2.78698E-03 

2.80556E-03 

2.82414E-03 

2.84272E-03 

2.86130E-03 

2.87988E-03 

2.89846E-03 

2.91704E-03 

2.93562E-03 

2.95420E-03 

2.97278E-03 

2.99136E-03 

3.00994E-03 

3.02852E-03 

3.04710E-03 

3.06568E-03 

3.08426E-03 

3.10284E-03 

3.12142E-03 

3.14000E-03 

3.15858E-03 

3.17716E-03 

3.19574E-03 

3.21432E-03 

3,23290E-03 

3.25148E-03 

3.27006E-03 

3.28864E-03 

3.30721E-03 

3.32579E-03 

3.34437E-03 

3.36295E-03 

3.38153E-03 

3,40011E-03 

3.41869E-03 

3.43727E-03 

3.45585E-03 

3,47443E-03 

3,49301E-03 

3,51159E-03 

3.53017E-03 

3.54875E-03 

3.56733E-03  3.58591E-03 

3.71597E-03 

Radial  distribution,  f(p,z). 

3.60449E-03 

cm'^ 

3.62307E-03 

3.64165E-03 

3.66023E-03 

3.67881E-03 

3.69739E-03 

2.41382E+06 

2.40689E+06 

2.38624E+06 

2.35222E+06 

2.30544E+06 

2.24671E+06 

2.17705E+06 

2.09763E+06 

2.00973E+06 

1.91476E+06 

1.81415E+06 

1.70937E+06 

1.60188E+06 

1.49306E+06 

1.38425E+06 

1.27665E+06 

1.17136E+06 

1.06934E+06 

9.71393E+05 

8.78184E+05 

7.90219E+05 

7.07863E+05 

6.31341E+05 

5.60757E+05 

4.96100E+05 

4.37265E+05 

3,84066E+05 

3.36252E+05 

2.93522E+05 

2.55541E+05 

2.21951E+05 

1.92384E+05 

1.66472E+05 

l,43854E+05 

1.24181E+05 

1.07126E+05 

9.23815E+04 

7.96666E+04 

6.87240E+04 

5.93224E+04 

5.12552E+04 

4.43393E+04 

3.84139E+04 

3.33382E+04 

2.89903E+04 

2.52645E+04 

2.20699E+04 

1.93285E+04 

1.69736E+04 

1.49480E+04 

1.32033E+04 

1.16981E+04 

1.03971E+04 

9.27032E+03 

8.29237E+03 

7.44156E+03 

6.69948E+03 

6.05048E+03 

5.48126E+03 

4.98051E+03 

4.53862E+03 

4.14741E+03 

3.79995E+03 

3.49032E+03 

3.21352E+03 

2.96526E+03 

2.74190E+03 

2,54033E+03 

2.35791E+03 

2.19237E+03 

2.04175E+03 

1.90437E+03 

1.77879E+03 

1.66375E+03 

1.55816E+03 

1.46105E+03 

1.37161E+03 

1.28908E+03 

1.21282E+03 

1.14226E+03 

1.07687E+03 

1.01620E+03 

9.59841E+02 

9.07423E+02 

8.58615E+02 

8.13118E+02 

7.70663E+02 

7.31004E+02 

6.93920E+02 

6.59210E+02 

6.26691E+02 

5.96196E+02 

5.67573E+02 

5.40684E+02 

5.15402E+02 

4.91611E+02 

4.69204E+02 

4.48086E+02 

4.28165E+02 

4.09361E+02 

3.91597E+02 

3.74804E+02 

3.58917E+02 

3.43879E+02 

3.29633E+02 

3.16130E+02 

3.03323E+02 

2.91169E+02 

2.79627E+02 

2.68660E+02 

2.58233E+02 

2.48315E+02 

2.38876E+02 

2.29887E+02 

2.21323E+02 

2.13160E+02 

2.05375E+02 

1.97946E+02 

1.90855E+02 

1.84082E+02 

1.77611E+02 

1.71426E+02 

1.65510E+02 

1.59850E+02 

1.54434E+02 

1.49247E+02 

1.44278E+02 

1.39517E+02 

l,34953E+02 

1.30576E+02 

1.26377E+02 

1.22346E+02 

1.18477E+02 

1.14762E+02 

1.11192E+02 

1.07761E+02 

1.04463E+02 

1.01291E+02 

9.82393E+01 

9.53031E+01 

9.24767E+01 

8.97552E+01 

8.71340E+01 

8.46085E+01 

8.21747E+01 

7.98284E+01 

7.75660E+01 

7.53838E+01 

7.32784E+01 

7,12465E+01 

6.92852E+01 

6.73913E+01 

6.55622E+01 

6.37952E+01 

6.20878E+01 

6.04375E+01 

5.88421E+01 

5,72993E+01 

5.58071E+01 

5.43635E+01 

5,29665E+01 

5.16145E+01 

5.03056E+01 

4.90382E+01 

4.78107E+01 

4.66216E+01 

4.54695E+01 

4.43529E+01 

4.32706E+01 

4.22214E+01 

4.12039E+01 

4.02171E+01 

3.92598E+01 

3.83310E+01 

3.74297E+01 

3.65548E+01 

3.57055E+01 

3.48808E+01 

3.40799E+01 

3,33021E+01 

3.25464E+01 

3.18121E+01 

3.10985E+01 

3.04049E+01 

2.97307E+01 

2.90751E+01 

2.84375E+01 

2.78175E+01 

2.72143E+01 

2.66275E+01 

2.60565E+01 

2.55007E+01 

2.49598E+01 

2.44332E+01 

2.39205E+01 

2.34212E+01 

2.29350E+01 

2.24613E+01 

2.19999E+01 

2.15503E+01 

2.11121E+01 
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Table  8.  Output  file  MORAD160.050  from  program  MORAD. 


Program  MORAD 

, output  file  MORAD160.050 

0 depths 
201  radial 

distances 

TIN  = energy  of  incident  proton  beam.  MeV 

R6IN  = CSDA 

range  at  energy  TIN,  g/cm 

ZR  = depth 

, in  units  of  RGIN 

RMAX  = largest  radial  distance 

SUM  = cumulative  integral 

of  Moliire  distribution  up  to 

RMAX 

TIN 

RGIN 

ZR 

RMAX 

SUM 

160.000000 

17.653381 

0.500000 

1.668557 

0.999355 

Radial  distances,  p,  cm 

O.OOOOOE+00 

8.34279E-03 

1.66856E-02 

2.S0284E-02 

3.33711E-02 

4.17139E-02 

5.00567E-02 

5.83995E-02 

6.67423E-02 

7.50851E-02 

8.34279E-02 

9.17707E-02 

1.00113E-01 

1.08456E-01 

1.16799E-01 

1.25142E-01 

1.33485E-01 

1.41827E-01 

1.50170E-01 

1.58513E-01 

1.66856E-01 

1.75199E-01 

1.83541E-01 

1.91884E-01 

2.00227E-01 

2.08570E-01 

2.16912E-01 

2.25255E-01 

2.33598E-01 

2.41941E-01 

2.50284E-01 

2.58626E-01 

2.66969E-01 

2.75312E-01 

2.83655E-01 

2.91998E-01 

3.00340E-01 

3.08683E-01 

3.17026E-01 

3.25369E-01 

3.33711E-01 

3.42054E-01 

3.50397E-01 

3.58740E-01 

3.67083E-01 

3.75425E-01 

3.83768E-01 

3.92111E-01 

4.00454E-01 

4.08797E-01 

4.17139E-01 

4.25482E-01 

4.33825E-01 

4.42168E-01 

4.50511E-01 

4.58853E-01 

4.67196E-01 

4.75539E-01 

4.83882E-01 

4.92224E-01 

5.00567E-01 

5.08910E-01 

5.17253E-01 

5.25596E-01 

5.33938E-01 

5.42281E-01 

5.50624E-01 

5.58967E-01 

5.67310E-01 

5.75652E-01 

5.83995E-01 

5.92338E-01 

6.00681E-01 

6.09023E-01 

6.17366E-01 

6.25709E-01 

6.34052E-01 

6.42395E-01 

6.50737E-01 

6.59080E-01 

6.67423E-01 

6.75766E-01 

6.84109E-01 

6.92451E-01 

7.00794E-01 

7.09137E-01 

7.17480E-01 

7.25822E-01 

7.34165E-01 

7.42508E-01 

7.50851E-01 

7.59194E-01 

7.67536E-01 

7.75879E-01 

7.84222E-01 

7.92565E-01 

8.00908E-01 

8.09250E-01 

8.17593E-01 

8.25936E-01 

8.34279E-01 

8.42622E-01 

8.50964E-01 

8.59307E-01 

8.67650E-01 

8.75993E-01 

8.84335E-01 

8.92678E-01 

9.01021E-01 

9.09364E-01 

9.17707E-01 

9.26049E-01 

9.34392E-01 

9.42735E-01 

9.51078E-01 

9.59421E-01 

9.67763E-01 

9.76106E-01 

9.84449E-01 

9.92792E-01 

1.00113E+00 

1.00948E+00 

1.01782E+00 

1.02616E+00 

1.03451E+00 

1.04285E+00 

1.05119E+00 

1.05953E+00 

1.06788E+00 

1.07622E+00 

1.08456E+00 

1.09291E+00 

1.10125E+00 

1.10959E+00 

1.11793E+00 

1.12628E+00 

1.13462E+00 

1.14296E+00 

1.15130E+00 

1.15965E+00 

1.16799E+00 

1.17633E+00 

1.18468E+00 

1.19302E+00 

1.20136E+00 

1.20970E+00 

1.21805E+00 

1.22639E+00 

1.23473E+00 

1.24308E+00 

1.25142E+00 

1.25976E+00 

1.26810E+00 

1.27645E+00 

1.28479E+00 

1.29313E+00 

1.30147E+00 

1.30982E+00 

1.31816E+00 

1.32650E+00 

1.33485E+00 

1.34319E+00 

1.35153E+00 

l,35987E+00 

1.36822E+00 

1.37656E+00 

1.38490E+00 

1.39325E+00 

1.40159E+00 

1.40993E+00 

1.41827E+00 

1.42662E+00 

1.43496E+00 

1.44330E+00 

1.45164E+00 

1.45999E+00 

1.46833E+00 

1.47667E+00 

1.48502E+00 

1.49336E+00 

1.50170E+00 

1.51004E+00 

1.51839E+00 

1.52673E+00 

1.53507E+00 

1.54342E+00 

1.55176E+00 

1.56010E+00 

1.56844E+00 

1.57679E+00 

1.58513E+00 

1.59347E+00 

1.60182E+00 

1.66856E+00 

1.61016E+00 

1.61850E+00 

1.62684E+00 

1.63519E+00 

1.64353E+00 

1.65187E+00 

1.66021E+00 

Radial  distribution,  f(p. 

z),  cm*^ 

1.18020E+01 

1.17695E+01 

1.16725E+01 

1.15127E+01 

1.12927E+01 

1.10163E+01 

1.06878E+01 

1.03127E+01 

9.89674E+00 

9.44622E+00 

8,96775E+00 

8.46803E+00 

7.95373E+00 

7.43134E+00 

6.90705E+00 

6.38660E+00 

5.87524E+00 

5.37761E+00 

4.89772E+00 

4.43889E+00 

4.00380E+00 

3.59445E+00 

3.21220E+00 

2.85784E+00 

2.53162E+00 

2.23332E+00 

1.96229E+00 

1.71759E+00 

1.49796E+00 

1.30198E+00 

1.12804E+00 

9.74488E-01 

8.39604E-01 

7.21677E-01 

6.19036E-01 

5.30071E-01 

4.53259E-01 

3.87176E-01 

3.30508E-01 

2.82055E-01 

2.40732E-01 

2.05568E-01 

1.75700E-01 

1.50369E-01 

1.28907E-01 

1.10737E-01 

9.53589E-02 

8.23421E-02 

7.13191E-02 

6.19764E-02 

5.40481E-02 

4.73094E-02 

4.15705E-02 

3.66719E-02 

3.24797E-02 

2.88817E-02 

2.57839E-02 

2.31077E-02 

2.07876E-02 

1.87685E-02 

1.70049E-02 

1.54582E-02 

1.40965E-02 

1.28930E-02 

1.18252E-02 

1.08742E-02 

1.00242E-02 

9.26176E-03 

8.57557E-03 

7.95604E-03 

7.39504E-03 

6.88559E-03 

6.42174E-03 

5.99836E-03 

5.61104E-03 

5.25594E-03 

4.92972E-03 

4.62947E-03 

4.35264E-03 

4.09696E-03 

3.86046E-03 

3.64136E-03 

3.43811E-03 

3.24929E-03 

3.07366E-03 

2.91010E-03 

2.75759E-03 

2.61523E-03 

2.48220E-03 

2.35774E-03 

2.24120E-03 

2.13196E-03 

2.02945E-03 

1.93319E-03 

1.84271E-03 

1.75758E-03 

1.67742E-03 

1.60188E-03 

1.53064E-03 

1.46340E-03 

1.39989E-03 

1.33985E-03 

1.28307E-03 

1.22931E-03 

1.17840E-03 

1.13014E-03 

1.08437E-03 

1.04093E-03 

9.99688E-04 

9.60500E-04 

9.23243E-04 

8.87805E-04 

8.54076E-04 

8.21959E-04 

7.91359E-04 

7.62191E-04 

7.34373E-04 

7.07830E-04 

6.82492E-04 

6.58293E-04 

6.35171E-04 

6.13068E-04 

5.91930E-04 

5.71707E-04 

5.52350E-04 

5.33816E-04 

5.16061E-04 

4.99047E-04 

4.82736E-04 

4.67093E-04 

4.52086E-04 

4.37683E-04 

4.23855E-04 

4.10574E-04 

3.97814E-04 

3.85552E-04 

3.73763E-04 

3.62425E-04 

3.51518E-04 

3.41022E-04 

3.30918E-04 

3.21189E-04 

3.11818E-04 

3.02790E-04 

2.94088E-04 

2.85700E-04 

2.77611E-04 

2.69809E-04 

2.62281E-04 

2.55015E-04 

2.48002E-04 

2.41230E-04 

2.34689E-04 

2.28370E-04 

2.22264E-04 

2.16362E-04 

2.10656E-04 

2.05139E-04 

1.99802E-04 

1.94638E-04 

1.89642E-04 

1.84806E-04 

1.80124E-04 

l,75590E-04 

1.71199E-04 

1.66945E-04 

1.62824E-04 

1.58829E-04 

1.54957E-04 

1.51204E-04 

1.47563E-04 

1.44033E-04 

1.40607E-04 

1.37284E-04 

1.34059E-04 

1.30928E-04 

1.27889E-04 

1.24938E-04 

1.22072E-04 

1.19289E-04 

1.16584E-04 

1.13956E-04 

1.11403E-04 

1.08920E-04 

1.06507E-04 

1.04161E-04 

1.01879E-04 

9.96594E-05 

9.75004E-05 

9.53999E-05 

9,33558E-05 

9.13665E-05 

8.94301E-05 

7.56552E-05 

8.75450E-05 

8.57096E-05 

8.39221E-05 

8.21813E-05 

8.04856E-05 

7.88335E-05 

7.72238E-05 
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Table  9.  Excerpt  from  output  file  PTRAD.160  from  program  PTRAD. 


Program  PTRAD,  output  file  PTRAD. 160 

9 depths 

201  radial  distances 

TIN  = energy  of  incident  proton  beam,  MeV 

BR  = depth,  in  units  of  CSDA  range  at  energy  TIN 

RTOP  = Maxinun  radial  distance 

TIN  ZR  RTOP 

160.000000  0.985000  2.648007 

Radial  distances,  p,  cm 

O.OOOOOE+00 

1.32A00E-02 

2.6A801E-02 

3.97201E-02 

5.29601E-02 

6.62002E-02 

7.9AA02E-02 

9.26803E-02 

1.05920E-01 

1.19160E-01 

1.32A00E-01 

1.A56A0E-01 

1.58880E-01 

1.72120E-01 

1.85361E-01 

1.98601E-01 

2,118A1E-01 

2.25081E-01 

2.38321E-01 

2.51561E-01 

2.6A801E-01 

2.780A1E-01 

2.91281E-01 

3.0A521E-01 

3.17761E-01 

3.31001E-01 

3.AA2A1E-01 

3.57A81E-01 

3.70721E-01 

3.83961E-01 

3.97201E-01 

A.lOAAlE-01 

A.23681E-01 

A.36921E-01 

A.50161E-01 

A.63A01E-01 

A.766A1E-01 

A.89881E-01 

5.03121E-01 

5.16361E-01 

5.29601E-01 

5.A28A1E-01 

5.56082E-01 

5.69322E-01 

5.82562E-01 

5.95802E-01 

6.090A2E-01 

6.22282E-01 

6.35522E-01 

6.A8762E-01 

6.62002E-01 

6.752A2E-01 

6.88A82E-01 

7.01722E-01 

7.1A962E-01 

7.28202E-01 

7.A1AA2E-01 

7.5A682E-01 

7.67922E-01 

7.81162E-01 

7.9AA02E-01 

8.076A2E-01 

8.20882E-01 

8.3A122E-01 

8.A7362E-01 

8.60602E-01 

8.738A2E-01 

8.87082E-01 

9.00322E-01 

9.13562E-01 

9.26803E-01 

9.A00A3E-01 

9.53283E-01 

9.66523E-01 

9.79763E-01 

9.93003E-01 

l,0062AE+00 

1.019A8E+00 

1.03272E+00 

1.0A596E+00 

1.05920E+00 

1.072AAE+00 

1.08568E+00 

1.09892E+00 

1.11216E+00 

1.125AOE+00 

1.1386AE+00 

1.15188E+00 

1.16512E+00 

1.17836E+00 

1.19160E+00 

1.20A8AE+00 

1.21808E+00 

1.23132E+00 

1.2AA56E+00 

1.25780E+00 

1.2710AE+00 

1.28A28E+00 

1,29752E+00 

1.31076E+00 

1.32A00E+00 

1.3372AE+00 

1.350A8E+00 

l,36372E+00 

1.37696E+00 

1.39020E+00 

1.A03AAE+00 

1.A1668E+00 

1.A2992E+00 

1.AA316E+00 

1.A56A0E+00 

1.A696AE+00 

1.A8288E+00 

1.A9612E+00 

1.50936E+00 

1.52260E+00 

1.5358AE+00 

l,5A908E+00 

1.56232E+00 

1.57556E+00 

1.58880E+00 

1.6020AE+00 

1.61528E+00 

1.62852E+00 

1.6A176E+00 

1.65500E+00 

1.6682AE+00 

1.681A8E+00 

1.69A72E+00 

1.70796E+00 

1.72120E+00 

1.73AAAE+00 

1.7A768E+00 

1.76092E+00 

1.77A16E+00 

1.787A0E+00 

1.8006AE+00 

1.81388E+00 

1,82712E+00 

1.8A036E+00 

1.85361E+00 

1.86685E+00 

1.88009E+00 

1.89333E+00 

1.90657E+00 

1.91981E+00 

1.93305E+00 

1.9A629E+00 

1.95953E+00 

1.97277E+00 

1.98601E+00 

1.99925E+00 

2.012A9E+00 

2.02573E+00 

2.03897E+00 

2.05221E+00 

2.065A5E+00 

2,07869E+00 

2,09193E+00 

2.10517E+00 

2.118A1E+00 

2.13165E+00 

2.1AA89E+00 

2.15813E+00 

2.17137E+00 

2.18A61E+00 

2.19785E+00 

2,21109E+00 

2.22A33E+00 

2.23757E+00 

2.25081E+00 

2.26A05E+00 

2.27729E+00 

2.29053E+00 

2.30377E+00 

2.31701E+00 

2.33025E+00 

2.3A3A9E+00 

2.35673E+00 

2.36997E+00 

2.38321E+00 

2.396A5E+00 

2.A0969E+00 

2.A2293E+00 

2.A3617E+00 

2.AA9A1E+00 

2.A6265E+00 

2.A7589E+00 

2.A8913E+00 

2.50237E+00 

2.51561E+00 

2.52885E+00 

2.5A209E+00 

2.55533E+00 

2.56857E+00 

2.58181E+00 

2.59505E+00 

2.60829E+00 

2.62153E+00 

2.63A77E+00 

2.6A801E+00 

Radial  distribution,  f(p,z) 

, cm'^ 

1.12791E+00 

1.12028E+00 

1.11265E+00 

1.10621E+00 

1.09929E+00 

1.09155E+00 

1.08288E+00 

1.0732AE+00 

1.06259E+00 

1.0509AE+00 

1.03826E+00 

1.02A55E+00 

1.00982E+00 

9.9A053E-01 

9.77252E-01 

9.59A16E-01 

9.A05AAE-01 

9.206A1E-01 

8.99821E-01 

8.78278E-01 

8.56168E-01 

8.33618E-01 

8,10733E-01 

7.87599E-01 

7.6A288E-01 

7.A0859E-01 

7.1736AE-01 

6.938A6E-01 

6.703A3E-01 

6.A6886E-01 

6.2350AE-01 

6.00220E-01 

5.77056E-01 

5.5A0A9E-01 

5.3126AE-01 

5.08762E-01 

A.86597E-01 

A.6A818E-01 

A,A3A67E-01 

A.22585E-01 

A.02205E-01 

3.82360E-01 

3.63078E-01 

3.AA385E-01 

3.2630AE-01 

3.08857E-01 

2.9206AE-01 

2.759A2E-01 

2.60501E-01 

2.A5723E-01 

2.31586E-01 

2.18070E-01 

2.05155E-01 

1.9282AE-01 

1.81061E-01 

1.69851E-01 

1.59180E-01 

1.A903AE-01 

1.39A02E-01 

1.30272E-01 

1.2163AE-01 

1.13A79E-01 

1.05796E-01 

9.85758E-02 

9.17983E-02 

8.5A38AE-02 

7.9A735E-02 

7.38813E-02 

6.86A17E-02 

6.37352E-02 

5.91A38E-02 

5,A8501E-02 

5.08381E-02 

A.7092AE-02 

A.35986E-02 

A.03A28E-02 

3.73121E-02 

3.AA9AAE-02 

3.18775E-02 

2.9A502E-02 

2.7200AE-02 

2.51173E-02 

2.31902E-02 

2.1A092E-02 

1.976A5E-02 

1.82A72E-02 

1.68A83E-02 

1.55598E-02 

1.A3736E-02 

1.32823E-02 

1.2278AE-02 

1.13552E-02 

1.05060E-02 

9.72AA3E-03 

9.00519E-03 

8.3A3A6E-03 

7.735A0E-03 

7.1767AE-03 

6.66382E-03 

6.1930AE-03 

5.76070E-03 

5.36350E-03 

A.99818E-03 

A.66172E-03 

A.35092E-03 

A.0630AE-03 

3.79517E-03 

3.5AA85E-03 

3.30926E-03 

3.08637E-03 

2.875A5E-03 

2.67619E-03 

2.A8815E-03 

2.31097E-03 

2.1AA31E-03 

1.98796E-03 

1.8A1A6E-03 

1.70A57E-03 

1.57698E-03 

1.A58A7E-03 

1.3A862E-03 

1.2A735E-03 

1.15A22E-03 

1.06913E-03 

9.91735E-0A 

9.21571E-0A 

8.58035E-0A 

8.0059AE-0A 

7,A8788E-0A 

7.01956E-OA 

6.59722E-0A 

6.21A60E-0A 

5.86822E-0A 

5.55318E-0A 

5.26A18E-0A 

A.9981AE-0A 

A.7A953E-0A 

A.515A9E-0A 

A.29126E-0A 

A.07368E-0A 

3.86018E-0A 

3.65165E-0A 

3.AA797E-0A 

3.25000E-0A 

3.05761E-0A 

2.87117E-0A 

2.69102E-0A 

2.517A8E-0A 

2.350A3E-0A 

2.19066E-0A 

2.03801E-0A 

1.89279E-0A 

1.75531E-0A 

1.62586E-0A 

1.50A73E-0A 

1.39131E-0A 

1.28633E-0A 

1.18833E-0A 

1.09717E-0A 

1.01316E-0A 

9.3A881E-05 

8.6307AE-05 

7.96359E-05 

7.35059E-05 

6.78660E-05 

6.26659E-05 

5.78979E-05 

5.35132E-05 

A.95051E-05 

A.58266E-05 

A.2A317E-05 

3.9355AE-05 

3.65129E-05 

3.39395E-05 

3.15913E-05 

2.9A258E-05 

2.7A783E-05 

2.57067E-05 

2.A0699E-05 

2.26036E-05 

2.12292E-05 

1.99830E-05 

1.88627E-05 

1.77918E-05 

1.68067E-05 

1.58690E-05 

1.A9780E-05 

1.A1330E-05 

1.33332E-05 

1.25A18E-05 

1.18305E-05 

1.11623E-05 

1.05011E-05 

9.88199E-06 

9.30A39E-06 

8.73271E-06 

8.20162E-06 

7.71050E-06 

7.22A3AE-06 

6.7A306E-06 

6.30065E-06 
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Table  10.  Ou^ut  file  FCIR.  1 from  program  FCIR. 


Program  FCIR,  outpxjt  file  FCIR.1 

Input  file  MORAD.160 

Input  file  PTRAD.160 

Input  file  PTPOL.160 

TIN  = energy  of  incident  proton  beam.  MeV 
RGIN  = CSOA  range  at  energy  TIN,  g/cnr 

RADIUS  = radius  of  circular  field,  cm 

LMAX  = number  of  radial  distances  from  center 
NCASE  = nunber  of  depths 

TIN  RGIN  RADIUS 

160.00000  17.65338  0.20000 

of  field 

LMAX 

1 

NCASE 

107 

Radial  distance  (cm)  from  center  of  field 

= 0.000 

Depths,  in 

units  of  RGIN 

0.010000 

0.020000 

0.030000 

0.040000 

0.050000 

0.060000 

0.070000 

0.080000 

0.090000 

0.100000 

0.110000 

0.120000 

0.130000 

0.140000 

0.150000 

0.160000 

0.170000 

0.180000 

0.190000 

0.200000 

0.210000 

0.220000 

0.230000 

0.240000 

0.250000 

0.260000 

0.270000 

0.280000 

0.290000 

0.300000 

0.310000 

0.320000 

0.330000 

0.340000 

0.350000 

0.360000 

0.370000 

0.380000 

0.390000 

0.400000 

O.AIOOOO 

0.420000 

0.430000 

0.440000 

0.450000 

0.460000 

0.470000 

0.480000 

0.A90000 

0.500000 

0.510000 

0.520000 

0.530000 

0.540000 

0.550000 

0.560000 

0.570000 

0.580000 

0.590000 

0.600000 

0.610000 

0.620000 

0.630000 

0.640000 

0.650000 

0.660000 

0.670000 

0.680000 

0.690000 

0.700000 

0.710000 

0.720000 

0.730000 

0.740000 

0.750000 

0.760000 

0.770000 

0.780000 

0.790000 

0.800000 

0.810000 

0.820000 

0.830000 

0.840000 

0.850000 

0.860000 

0.870000 

0.880000 

0.890000 

0.900000 

0.910000 

0.920000 

0.930000 

0.940000 

0.950000 

0.960000 

0.970000 

0.980000 

0.985000 

0.990000 

0.995000 

1.000000 

1.005000 

1.010000 

1.015000 

1.020000 

1.025000 

Reduction  factors 

1.000000 

1.000000 

1.000000 

1.000000 

1.000000 

1.000000 

1.000000 

1.000000 

1.000000 

1.000000 

1.000000 

1.000000 

1.000000 

0.999832 

0.999614 

0.999361 

0.999070 

0.998734 

0.998350 

0.997912 

0.997409 

0.996836 

0.996181 

0.995431 

0.994569 

0.993575 

0.992419 

0.991062 

0.989448 

0.987506 

0.985141 

0.982245 

0.978695 

0.974354 

0.969090 

0.962774 

0.955300 

0.946574 

0.936539 

0.925168 

0.912465 

0.898469 

0.883240 

0.866869 

0.849459-^ 

^ 0.831132 

0.812014 

0.792238 

0.771932 

0.751228 

0.730247 

0.709103 

0.687904 

0.666744 

0.645713 

0.624881 

0.604320 

0.584085 

0.564225 

0.544779 

0.525781 

0.507255 

0.489221 

0.471694 

0.454684 

0.438194 

0.422227 

0.406782 

0.391855 

0.377439 

0.363525 

0.350105 

0.337167 

0.324700 

0.312691 

0.301126 

0.289991 

0.279274 

0.268961 

0.259035 

0.249485 

0.240295 

0.231453 

0.222944 

0.214755 

0.206874 

0.199286 

0.191981 

0.184944 

0.178166 

0.171632 

0.165332 

0.159254 

0.153386 

0.147717 

0.142232 

0.136919 

0.131760 

0.129551 

0.128709 

0.127830 

0.128633 

0.130434 

0.133959 

0.133396 

0.139213 

0.139700 

Absorbed  dose,  NeV/g,  from  a beam  of  1 proton/cm^ 

6.249280 

6.253005 

6.256818 

6.260764 

6.264886 

6.269230 

6.273839 

6.278757 

6.284030 

6.289700 

6.295804 

6.302346 

6.309319 

6.315659 

6.322102 

6.328735 

6.335532 

6.342449 

6.349453 

6.356497 

6.363521 

6.370490 

6.377352 

6.384037 

6.390450 

6.396480 

6.401947 

6.406600 

6.410099 

6.411972 

6.411624 

6.408335 

6.401297 

6.389614 

6.372371 

6.348702 

6.317856 

6.279141 

6.232118 

6.176513 

6.112279 

6.039575 

5.958720 

5.870231 

5.774710 

5.672882 

5.565526 

5.453476 

5.337544 

5.218554 

5.097287 

4.974484 

4.850840 

4.726976 

4.603480 

4.480829 

4.359514 

4.239927 

4.122427 

4.007289 

3.894760 

3.785048 

3.678329 

3.574723 

3.474328 

3.377211 

3.283452 

3.193109 

3.106209 

3.022755 

2.942734 

2.866135 

2.792947 

2.723197 

2.656898 

2.594047 

2.534628 

2.478673 

2.426200 

2.377294 

2.332048 

2.290564 

2.252980 

2.219541 

2.190535 

2.166381 

2.147519 

2.134632 

2.128588 

2.131040 

2.143681 

2.168993 

2.211883 

2.279476 

2.385033 

2.558756 

2.877361 

3.413780 

3.686637 

3.824325 

3.620271 

3.064799 

2.216462 

1.356335 

0.655387 

0.267163 

0.082841 
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Table  11.  Output  file  FCIR.2  from  program  FCIR. 


Program  FCIR,  output  file  FCIR.2 

Input  file  MORAD160.050 

Input  file  PTPOL.160 

TIN  = energy  of  incident  proton  beam.  MeV 
RGIN  = CSDA  range  at  energy  TIN,  g/air 

RADIUS  = radius  of  circular  field,  cm 

LHAX  = number  of  radial  distances  from  center 
NCASE  = number  of  depths 

TIN  RGIN  RADIUS 

160.00000  17.65338  0.20000 

of  field 

LMAX 

101 

NCASE 

1 

Radial  distances  from  center  of  field  (cm) 

0.000000 

0.005000 

0.010000 

0.015000 

0.020000 

0.025000 

0.030000 

0.035000 

O.OAOOOO 

0.045000 

0.050000 

0.055000 

0.060000 

0.065000 

0.070000 

0.075000 

0.080000 

0.085000 

0.090000 

0.095000 

0.100000 

0.105000 

0.110000 

0.115000 

0.120000 

0.125000 

0.130000 

0.135000 

0.140000 

0.145000 

0.150000 

0.155000 

0.160000 

0.165000 

0.170000 

0.175000 

0.180000 

0.185000 

0.190000 

0.195000 

0.200000 

0.205000 

0.210000 

0.215000 

0.220000 

0.225000 

0.230000 

0.235000 

0.240000 

0.245000 

0.250000 

0.255000 

0.260000 

0.265000 

0.270000 

0.275000 

0.280000 

0.285000 

0.290000 

0.295000 

0.300000 

0.305000 

0.310000 

0.315000 

0.320000 

0.325000 

0.330000 

0.335000 

0.340000 

0.345000 

0.350000 

0.355000 

0.360000 

0.365000 

0.370000 

0.375000 

0.380000 

0.385000 

0.390000 

0.395000 

0.400000 

0.405000 

0.410000 

0.415000 

0.420000 

0.425000 

0.430000 

0.435000 

0.440000 

0.445000 

0.450000 

0.455000 

0.460000 

0.465000 

0.470000 

0.475000 

0.480000 

0.485000 

0.490000 

0.495000 

0.500000 

Depth  as 

fraction  of  RGIN 

= 0.500 

Reduction 

factors 

0.751228 

0.750933 

0.750050 

0.748579 

0.746523 

0.743882 

0.740661 

0.736862 

0.732492 

0.727554 

0.722055 

0.716003 

0.709404 

0.702268 

0.694606 

0.686427 

0.677744 

0.668571 

0.658922 

0.648812 

0.638258 

0.627279 

0.615892 

0.604120 

0.591983 

0.579504 

0.566707 

0.553617 

0.540258 

0.526659 

0.512845 

0.498846 

0.484689 

0.470403 

0.456019 

0.441564 

0.427070 

0.412564 

0.398075 

0.383632 

0.369254 

0.354992 

0.340867 

0.326899 

0.313114 

0.299535 

0.286184 

0.273083 

0.260252 

0.247708 

0.235468 

0.223548 

0.211960 

0.200716 

0.189826 

0.179298 

0.169139 

0.159354 

0.149947 

0.140918 

0.132269 

0.123998 

0.116102 

0.108578 

0.101421 

0.094625 

0.088182 

0.082085 

0.076324 

0.070891 

0.065774 

0.060965 

0.056450 

0.052219 

0.048261 

0.044563 

0.041113 

0.037899 

0.034911 

0.032134 

0.029560 

0.027175 

0.024968 

0.022930 

0.021049 

0.019315 

0.017718 

0.016250 

0.014901 

0.013663 

0.012528 

0.011488 

0.010535 

0.009664 

0.008868 

0.008140 

0.007476 

0.006869 

0.006316 

0.005811 

0.005350 

Absorbed  dose 

, MeV/g,  from 

a beam  of  1 

2 

proton/ cm 

5.218554 

5.216508 

5.210374 

5.200157 

5.185868 

5.167525 

5.145148 

5.118762 

5.088400 

5.054099 

5.015901 

4.973855 

4.928016 

4.878447 

4.825216 

4.768401 

4.708086 

4.644363 

4.577331 

4.507100 

4.433787 

4.357516 

4.278420 

4.196642 

4.112331 

4.025644 

3.936746 

3.845810 

3.753013 

3.658541 

3.562582 

3.465332 

3.366988 

3.267751 

3.167826 

3.067416 

2.966725 

2.865957 

2.765311 

2.664979 

2.565095 

2.466025 

2.367898 

2.270869 

2.175107 

2.080778 

1.988036 

1.897029 

1.807892 

1.720754 

1.635728 

1.552919 

1.472419 

1.394311 

1.318661 

1.245529 

1.174960 

1.106987 

1.041636 

0.978917 

0.918832 

0.861374 

0.806526 

0.754259 

0.704541 

0.657329 

0.612572 

0.570216 

0.530200 

0.492457 

0.456916 

0.423503 

0.392142 

0.362752 

0.335253 

0.309563 

0.285598 

0.263275 

0.242513 

0.223229 

0.205342 

0.188773 

0.173446 

0.159285 

0.146218 

0.134173 

0.123084 

0.112885 

0.103515 

0.094915 

0.087028 

0.079802 

0.073187 

0.067135 

0.061602 

0.056547 

0.051931 

0.047717 

0.043872 

0.040366 

0.037168 
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Table  12.  Output  file  FREC.l  from  program  FREC. 


Program  FREC 

, output  file  FREC.l 

Input  file  = 

MORAD.160 

Input  file  = 

PTRAD.160 

Input  file  = 

PTPOL.160 

TIN 

= energy  of  incident  proton  beam.  MeV 

RGIN 

= CSDA  range  at  energy  TIN,  g/car 

2A 

- side  of  rectangular  field,  cm 

2B 

= side  of  rectangular  field,  cm 

ALPHA 

= angle  (deg)  defined  below 

LMAX 

= number  of  radial  distances  from  center 

of  field 

NCASE 

= number  of  depths 

Rectangular  field  extends  from  -A  to  A in  x,  and  from 

-B  to  B in  y. 

Line  along  which 

reduction 

factor  is 

calculated  starts  at  the  origin  and  makes  angle  ALPHA  (deg)  with  respect  to  x-axis. 

TIN 

RGIN  A 

B 

ALPHA 

LMAX 

NCASE 

160.00000  17.65338  0.05000 

0.50000 

0.00000 

1 

107 

Radial  distance  (cm)  from  center  of  field  = 0.000 
Depths,  in  units  of  RGIN 

0.010000 

0.020000  0.030000  0.040000 

0.050000 

0.060000 

0.070000 

0.080000 

0.090000 

0.100000  0.110000  0.120000 

0.130000 

0.140000 

0.150000 

0.160000 

0.170000 

0.180000  0.190000  0.200000 

0.210000 

0.220000 

0.230000 

0.240000 

0.250000 

0.260000  0.270000  0.280000 

0.290000 

0.300000 

0.310000 

0.320000 

0.330000 

0.340000  0.350000  0.360000 

0.370000 

0.380000 

0.390000 

0.400000 

0.410000 

0.420000  0.430000  0.440000 

0.450000 

0.460000 

0.470000 

0.480000 

0.490000 

0.500000  0.510000  0.520000 

0.530000 

0.540000 

0.550000 

0.560000 

0.570000 

0.580000  0.590000  0.600000 

0.610000 

0.620000 

0.630000 

0.640000 

0.650000 

0.660000  0.670000  0.680000 

0.690000 

0.700000 

0.710000 

0.720000 

0.730000 

0.740000  0.750000  0.760000 

0.770000 

0.780000 

0.790000 

0.800000 

0.810000 

0.820000  0.830000  0.840000 

0.850000 

0.860000 

0.870000 

0.880000 

0.890000 

0.900000  0.910000  0.920000 

0.930000 

0.940000 

0.950000 

0.960000 

0.970000 

0.980000  0.985000  0.990000 

0.995000 

1.000000 

1.005000 

1.010000 

1.015000 

1.020000  1.025000 

Reduction 

factors 

1.000000 

1.000000 

1.000000 

1.000000 

1.000000 

0.999921 

0.999655 

0 . 999228 

0.998628 

0.997801 

0.996676 

0.995124 

0.992895 

0.989566 

0.984554 

0 . 977221 

0.967035 

0.953648 

0.937119 

0.917651 

0.895618 

0.871597 

0.846089 

0.819596 

0.792565 

0.765385 

0.738368 

0.711759 

0.685743 

0.660457 

0.635995 

0.612418 

0.589760 

0.568035 

0.547274 

0.527396 

0.508415 

0.490302 

0.473025 

0.456551 

0.440844 

0.425868 

0.411586 

0.397963 

0.384966 

0.372560 

0.360713 

0.349395 

0.338577 

0.328230 

0.318328 

0.308880 

0.299792 

0.291077 

0.282713 

0.274680 

0.266957 

0.259528 

0.252373 

0.245476 

0.238822 

0.232395 

0.226183 

0.220170 

0.214347 

0.208700 

0.203220 

0.197897 

0 . 192722 

0.187686 

0.182783 

0.178006 

0.173349 

0.168805 

0.164371 

0.160041 

0.155811 

0.151678 

0.147638 

0.143689 

0.139828 

0.136051 

0.132358 

0.128745 

0.125210 

0.121752 

0.118369 

0.115058 

0.111818 

0.108648 

0.105544 

0.102504 

0.099528 

0.096611 

0.093752 

0.090946 

0.088189 

0.085475 

0.084419 

0.083850 

0.083572 

0.083792 

0.084856 

0.086104 

0.086641 

0.089510 

0.089876 

Absorbed 

dose,  MeV/g,  from 

a beam  of  1 

proton/ cm^ 

6.249280 

6.253005 

6.256818 

6.260764 

6.264886 

6.268737 

6.271671 

6.273911 

6.275405 

6.275869 

6.274878 

6.271616 

6.264490 

6.250812 

6.226854 

6.188526 

6.132385 

6.056129 

5.960025 

5.845253 

5.714090 

5.570127 

5.416494 

5.256345 

5.092504 

4.927431 

4.763102 

4.601079 

4.442557 

4.288415 

4.139266 

3.995518 

3.857411 

3.725052 

3.598672 

3.477744 

3.362391 

3.252439 

3.147705 

3.047982 

2.953055 

2.862714 

2.776737 

2.694915 

2.617036 

2.542902 

2.472319 

2.405110 

2.341098 

2.280117 

2.222005 

2.166848 

2.114020 

2.063629 

2.015546 

1.969641 

1.925807 

1.883934 

1.843924 

1.805674 

1.769091 

1.734095 

1.700610 

1.668556 

1.637865 

1.608474 

1.580340 

1.553424 

1.527691 

1.503105 

1.479631 

1.457249 

1.435945 

1.415736 

1.396639 

1.378669 

1.361842 

1.346202 

1.331793 

1.318707 

1.307033 

1.296881 

1.288378 

1.281729 

1.277162 

1.274987 

1.275548 

1.279330 

1.286955 

1.299533 

1.318240 

1.344756 

1.382343 

1.435740 

1.513715 

1.636112 

1.853288 

2.214559 

2.402326 

2.491435 

2.366847 

1.996419 

1.441960 

0.871803 

0.425676 

0.171778 

0.053296 
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Table  13.  Output  file  FREC.2A  from  program  FREC. 


Program  FREC,  output  file  FREC.2A 
Input  file  = HORAD160.050 
Input  file  = PTPOL.160 

TIN  = energy  of  incident  proton  beam.  MeV 
RGIN  = CSDA  range  at  energy  TIN,  g/cm^ 

2A  = side  of  rectangular  field,  cm 
2B  = side  of  rectangular  field,  cm 
ALPHA  = angle  (deg)  defined  below 
LMAX  = number  of  radial  distances  from  center  of  field 
NCASE  = number  of  depths 


Rectangular  field  extends  from  -A  to  . 

A in  X,  and  from 

-B  to  B in 

y.  Line  along  which 

reduction 

factor  is 

calculated  starts  at  the  origin  and  makes  angle  ALPHA 

(deg)  with 

respect  to  x-axis. 

TIN 

RGIN 

A 

B 

ALPHA 

LHAX 

NCASE 

160 

.00000 

17.65338 

0.05000 

0.50000 

0.00000 

101 

1 

Radial  distances  from  center 

of  field 

(cm) 

0.00000 

0.00500 

0.01000 

0.01500 

0.02000 

0.02500 

0.03000 

0.03500 

0.04000 

0.04500 

0.05000 

0.05500 

0.06000 

0.06500 

0.07000 

0.07500 

0.08000 

0.08500 

0.09000 

0.09500 

0.10000 

0.10500 

0.11000 

0.11500 

0.12000 

0.12500 

0.13000 

0.13500 

0.14000 

0.14500 

0.15000 

0.15500 

0.16000 

0.16500 

0.17000 

0.17500 

0.18000 

0.18500 

0.19000 

0.19500 

0.20000 

0.20500 

0.21000 

0.21500 

0.22000 

0.22500 

0.23000 

0.23500 

0.24000 

0.24500 

0.25000 

0.25500 

0.26000 

0.26500 

0.27000 

0.27500 

0.28000 

0.28500 

0.29000 

0.29500 

0.30000 

0.30500 

0.31000 

0.31500 

0.32000 

0.32500 

0.33000 

0.33500 

0.34000 

0.34500 

0.35000 

0.35500 

0.36000 

0.36500 

0.37000 

0.37500 

0.38000 

0.38500 

0.39000 

0.39500 

0.40000 

0.40500 

0.41000 

0.41500 

0.42000 

0.42500 

0.43000 

0.43500 

0.44000 

0.44500 

0.45000 

0.45500 

0.46000 

0.46500 

0.47000 

0.47500 

0.48000 

0.48500 

0.49000 

0.49500 

0.50000 

Depth  as 

fraction  of  RGIN 

= 0.500 

Reduction 

factors 

0.328230 

0.328015 

0.327025 

0.325536 

0.323571 

0.320883 

0.317758 

0.314051 

0.309860 

0.305112 

0.299950 

0.294321 

0.288282 

0.281896 

0.275192 

0.268123 

0.260767 

0.253223 

0.245407 

0.237490 

0.229369 

0.221173 

0.212943 

0.204630 

0.196290 

0.187997 

0.179778 

0.171636 

0.163557 

0.155617 

0.147874 

0.140275 

0.132826 

0.125571 

0.118607 

0.111781 

0.105230 

0.098903 

0.092867 

0.087030 

0.081471 

0.076161 

0.071123 

0.066311 

0.061760 

0.057458 

0.053395 

0.049559 

0.045951 

0.042581 

0.039398 

0.036431 

0.033654 

0.031077 

0.028664 

0.026427 

0.024348 

0.022426 

0.020641 

0.018993 

0.017471 

0.016070 

0.014775 

0.013586 

0.012493 

0.011488 

0.010565 

0.009719 

0.008946 

0.008235 

0.007585 

0.006990 

0.006448 

0.005949 

0.005494 

0.005078 

0.004699 

0.004351 

0.004034 

0.003743 

0.003478 

0.003235 

0.003012 

0.002808 

0.002622 

0.002451 

0.002293 

0.002149 

0.002017 

0.001894 

0.001782 

0.001678 

0.001583 

0.001494 

0.001412 

0.001336 

0.001266 

0.001200 

0.001140 

0.001083 

0.001031 

Absorbed  dose 

, MeV/g,  from 

a beam  of  1 

proton/cm^ 

2.280117 

2.278619 

2.271747 

2.261404 

2.247748 

2.229079 

2.207372 

2.181619 

2.152504 

2.119521 

2.083663 

2.044559 

2.002610 

1.958250 

1.911678 

1.862573 

1.811469 

1.759064 

1.704766 

1.649769 

1.593355 

1.536425 

1.479250 

1.421505 

1.363570 

1.305958 

1.248861 

1.192306 

1.136179 

1.081023 

1.027233 

0.974447 

0.922699 

0.872304 

0.823928 

0.776506 

0.731002 

0.687053 

0.645118 

0.604571 

0.565955 

0.529069 

0.494071 

0.460646 

0.429029 

0.399142 

0.370917 

0.344270 

0.319210 

0.295801 

0.273686 

0.253077 

0.233786 

0.215884 

0.199118 

0.183578 

0.169135 

0.155790 

0.143384 

0.131938 

0.121367 

0.111632 

0.102640 

0.094375 

0.086782 

0.079804 

0.073395 

0.067516 

0.062145 

0.057206 

0.052693 

0.048556 

0.044789 

0.041326 

0.038168 

0.035277 

0.032645 

0.030227 

0.028022 

0.026003 

0.024163 

0.022471 

0.020926 

0.019510 

0.018216 

0.017023 

0.015932 

0.014928 

0.014009 

0.013159 

0.012377 

0.011656 

0.010993 

0.010377 

0.009809 

0.009281 

0.008794 

0.008340 

0.007918 

0.007525 

0.007160 
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Table  14.  Output  file  FREC.2B  from  program  FREC. 


Program  FREC,  output  file  FREC.2B 
Input  file  = MORAD160.050 
Input  file  = PTPOL.160 

TIN  = energy  of  incident  proton  beam,  MeV 
RGIN  - CSDA  range  at  energy  TIN,  g/cin2 
2A  - side  of  rectangular  field,  cm 
2B  - side  of  rectangular  field,  cm 
ALPHA  = angle  (deg)  defined  below 
LMAX  = number  of  radial  distances  from  center  of  field 
NCASE  s number  of  depths 


Rectangular  field  extends  from  -A  to  A in  x,  and  from  -B  to  B in 
calculated  starts  at  the  origin  and  makes  angle  ALPHA  (deg)  with 

y.  Line  along  which 
respect  to  x-axis. 

reduction 

factor  is 

160 

TIN  RGIN 

.00000  17.65338 

A 

0.05000 

B 

0.50000 

ALPHA 

90.00000 

LHAX 

101 

NCASE 

1 

Radial  distances  from  center 

of  field 

(cm) 

0.00000 

0.01000 

0.02000 

0.03000 

0.04000 

0.05000 

0.06000 

0.07000 

0.08000 

0.09000 

0.10000 

0.11000 

0.12000 

0.13000 

0.14000 

0.15000 

0.16000 

0.17000 

0.18000 

0.19000 

0.20000 

0.21000 

0.22000 

0.23000 

0.24000 

0.25000 

0.26000 

0.27000 

0.28000 

0.29000 

0.30000 

0.31000 

0.32000 

0.33000 

0.34000 

0.35000 

0.36000 

0.37000 

0.38000 

0.39000 

0.40000 

0.41000 

0.42000 

0.43000 

0.44000 

0.45000 

0.46000 

0.47000 

0.48000 

0.49000 

0.50000 

0.51000 

0.52000 

0.53000 

0.54000 

0.55000 

0.56000 

0.57000 

0.58000 

0.59000 

0.60000 

0.61000 

0.62000 

0.63000 

0.64000 

0.65000 

0.66000 

0.67000 

0.68000 

0.69000 

0.70000 

0.71000 

0.72000 

0.73000 

0.74000 

0.75000 

0.76000 

0.77000 

0.78000 

0.79000 

0.80000 

0.81000 

0.82000 

0.83000 

0.84000 

0.85000 

0.86000 

0.87000 

0.88000 

0.89000 

0.90000 

0.91000 

0.92000 

0.93000 

0.94000 

0.95000 

0.96000 

0.97000 

0.98000 

0.99000 

1.00000 

Depth  as 

fraction  of  RGIN 

= 0.500 

Reduction 

factors 

0.328230 

0 . 328229 

0.328224 

0.328214 

0.328201 

0.328181 

0.328157 

0.328126 

0.328088 

0.328042 

0.327983 

0.327914 

0.327828 

0.327724 

0.327600 

0.327447 

0.327261 

0.327035 

0.326762 

0.326433 

0.326030 

0.325545 

0.324967 

0.324250 

0.323385 

0.322377 

0.321203 

0.319782 

0.318102 

0.316329 

0.313851 

0.311203 

0.308151 

0.304881 

0.300708 

0.296197 

0.291230 

0.285829 

0.279344 

0.272575 

0.265176 

0.257154 

0.249037 

0.239656 

0.229853 

0.219617 

0.209126 

0.198160 

0.187044 

0.175709 

0.164298 

0.152863 

0.141542 

0.130401 

0.119521 

0.108981 

0.098760 

0.089171 

0.080033 

0.071431 

0.063425 

0.056034 

0.049222 

0.043032 

0.037428 

0.032406 

0.027856 

0.023889 

0.020404 

0.017362 

0.014718 

0.012431 

0.010471 

0.008800 

0.007382 

0.006186 

0.005180 

0.004337 

0.003634 

0.003049 

0.002563 

0.002160 

0.001826 

0.001550 

0.001321 

0.001132 

0.000975 

0.000844 

0.000735 

0.000644 

0.000567 

0.000502 

0.000447 

0.000400 

0.000360 

0.000325 

0.000295 

0.000269 

0.000246 

0.000226 

0.000208 

Absorbed  dose 

, MeV/g,  from 

a beam  of  1 

proton/ cm2 

2.280117 

2.280109 

2.280071 

2.280007 

2.279913 

2.279778 

2.279610 

2.279393 

2.279131 

2.278807 

2.278399 

2.277919 

2.277326 

2.276600 

2.275738 

2.274677 

2.273386 

2.271815 

2.269919 

2.267633 

2.264832 

2.261464 

2.257448 

2.252464 

2.246455 

2.239459 

2.231304 

2.221432 

2.209758 

2.197440 

2.180229 

2.161835 

2.140632 

2.117918 

2.088932 

2.057592 

2.023090 

1.985567 

1.940517 

1.893495 

1.842101 

1.786369 

1.729988 

1.664815 

1.596716 

1.525617 

1.452736 

1.376558 

1.299340 

1.220600 

1.141332 

1.061892 

0.983249 

0.905858 

0.830274 

0.757059 

0.686054 

0.619443 

0.555968 

0.496207 

0.440593 

0.339252 

0.341930 

0.298928 

0.260002 

0.225114 

0.193507 

0.165949 

0.141741 

0.120609 

0 . 102243 

0.086355 

0.072739 

0.061130 

0.051284 

0.042971 

0.035982 

0.030129 

0.025243 

0.021177 

0.017801 

0.015002 

0.012685 

0.010767 

0.009179 

0.007864 

0.006772 

0.005864 

0.005107 

0.004473 

0.003940 

0.003489 

0.003106 

0.002781 

0.002501 

0.002261 

0.002053 

0.001871 

0.001712 

0.001572 

0.001448 
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Table  15.  Output  file  FRECl.l  from  program  FRECl. 


Program  FREC1,  output  file  FRECl.l 
Input  file  = HORAD50.160 
Input  file  = PTPOL.160 

Rectangular  field  extends  from  -A  to  A in  x,  and  from  -B  to  B in  y. 


TIN  = energy  of  incident  proton  beam.  NeV 
RGIN  = CSDA  range  at  energy  TIN,  g/avr 
2A  = side  of  rectangular  field,  cm 
2B  = side  of  rectangular  field,  cm 
X = x-distance  from  center  of  field,  cm 

Y = y-distance  from  center  of  field,  cm 

LMAX  = nunber  of  (X,Y)  values 

NCASE  = number  of  depths 

REDUC  = reduction  factor 

DRAD  = absorbed  dose,  MeV/g 


TIN 

160.00000 


RGIN 

17.65338 


A 

0.50000 


B 

0.20000 


LMAX 

9 


NCASE 

1 


Depths  as  fraction  of  RGIN  = 0.500 


X 

y 

O.AOOOO 

0.00000 

O.AOOOO 

0.20000 

O.AOOOO 

O.AOOOO 

0.50000 

0.00000 

0.50000 

0.20000 

0.50000 

O.AOOOO 

0.60000 

0.00000 

0.60000 

0.20000 

0.60000 

O.AOOOO 

REDUC 

DRAD 

0.72156 

5.012A6 

0.39753 

2.76153 

0.03857 

0.26795 

0.AA855 

3.11596 

0.2A801 

1.72282 

0.02505 

0.17A02 

0.17550 

1.21915 

0.09838 

0.68339 

0.01153 

0.08007 
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Figure  Captions 

Fig.  1.  Diagram  illustrating  the  calculation  of  the  reduction  factor.  The  solid  line  is  the 

boundary  of  the  irradiated  field.  The  fraction  of  the  circle  of  radius  p around  point  P 
that  lies  inside  the  irradiated  field  is  equal  to  the  sum  of  the  arc  lengths  AB  and  CD, 
divided  by  I'kq. 

Fig.  2.  Diagram  illustrating  the  calculation  of  the  reduction  factor  for  a rectangular  field.  The 
irradiated  field  is  the  rectangle  ABCD.  The  absoibed  dose  at  point  inside  the 
rectangle  can  be  calculated  as  the  sum  of  the  contributions  from  the  rectangles  P^AEF, 
PiFBG,  PjGCH  and  P^HDE.  The  absorbed  dose  at  point  P2  outside  the  rectangle  can  be 
calculated  as  the  sum  of  the  contributions  from  the  rectangles  P2KAI  and  P2JDK,  minus 
the  contributions  from  the  rectangles  P2LBI  and  P2LCJ.  (adapt^  from  Meredith  and 
Neary,  1944). 

Fig.  3.  Energy  deposition  distribution  dD/dz  for  proton  beams  with  energies  between  250  MeV 
and  50  MeV.  (From  Berger,  1993b). 

Fig.  4.  Scaled  radial  distributions  calculated  with  PTRAN.  The  solid  curves  are  for  a beam 
energy  of  160  MeV,  and  the  dotted  curves  for  70  MeV. 

a.  Results  at  shallow  and  intermediate  depths. 

b.  Results  at  great  depths,  at  and  beyond  the  Bragg  peak. 

Fig.  5.  Effect  of  energy-loss  straggling  on  the  radial  distributions  calculated  with  PTRAN.  The 
solid  curves  represent  results  that  take  into  account  energy-loss  straggling,  whereas  the 
points  (o)  represent  results  calculated  in  the  continuous-slowing-down  approximation. 

Fig.  6.  Percentage  amount  by  which  the  average  path  length  s^^  exceeds  the  penetration  depth, 
when  a proton  has  reach  depth  z.  The  squares  represent  Monte  Carlo  results  from 
PTRAN,  and  the  curve  is  from  a least-squares  cubic-spline  fit. 

Fig.  7.  Radial  dose  distribution,  27rp  f(p,z),  as  a function  of  the  radial  distance  from  the  axis  of  a 
narrow  pencil  beam.  The  histograms  were  calculated  with  PTRAN,  and  the  curves  were 
obtained  from  the  theory  of  Moli^re  (1955).  The  distributions  shown  are  normalized  to 
unity,  when  integrated  with  respect  to  p. 

a.  Beam  energy  T^  = 160  MeV 

b.  Beam  energy  Tq  = 70  MeV 

Fig.  8.  Moli^re  distribution  fM(^9)  as  function  of  function  of  the  scaled  angular  variable  t?,  for 
selected  values  of  the  parameter  B.  The  solid  curves  represent  the  Moli^re  theory,  and 
the  dotted  curves  represent  the  Gaussian  approximation. 

Fig.  9.  Reductions  factors  for  a circular  field,  at  depth  z/r^  = 0.1,  0.5  and  0.9,  for  a beam 


energy 

of  160  MeV. 

a. 

Field  radius  8 mm 

b. 

Field  radius  4 mm 

c. 

Field  radius  2 mm 

d. 

Field  radius  1 mm 
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Fig.  10.  Reduction  factors  along  the  central  axis,  for  a beam  energy  of  160  MeV  and  a circular 
field  with  radii  of  1,  2,  4 and  8 mm. 

Fig.  11.  Absorbed-dose  values  along  the  central  axis,  for  a beam  energy  of  160  MeV  and  a 
circular  field  with  radii  of  1 , 2 and  4 mm.  The  dotted  curve  pertains  to  a laterally 
unbounded  broad  beam. 

Fig.  12.  Absorbed-dose  values  along  lines  that  are  parallel  to  the  z axis  and  at  distances  of  0,2, 
3.98,  4.02,  6 and  8 mm  from  the  z axis,  for  a beam  energy  of  160  MeV. 

a.  Field  radius  4 mm 

b.  Field  radius  2 mm 

c.  Field  radius  1 mm 

Fig.  13.  Average  absorbed  dose  as  a function  of  depth,  for  circular  fields  with  radii  of  8,  4,  2 and 
1 mm.  The  absorbed  dose  is  averaged  over  thin  disks  with  radii  equal  to  the  field  radii. 

Fig.  14.  Reduction  factors  for  a circular  filed,  with  a pencil  beam  density  n(r)  = exp(-Qr^),  where 
r is  the  radial  distance  from  the  center  of  the  field.  The  results  are  for  a 160-MeV  beam 
and  a field  radius  of  4 mm. 

a.  Depth  z = 0. 1 r^ 

b.  Depth  z = 0.5  r^ 

c.  Depth  z = 0.9  r^ 

Fig.  15.  Comparison  of  reduction  factors  for  circular  and  square  fields  with  the  same  area 

(64  mm^),  for  a 160-MeV  beam.  The  curves  with  the  label  ”0  deg”  pertain  to  reduction 
factors  from  a square  field  along  a line  that  starts  a the  center  is  parallel  to  an  edge  of  the 
field.  The  label  ”45  deg"  pertains  to  a line  starting  at  the  center  that  passes  through  a 
comer  of  the  field. 

a.  Depth  z = 0. 1 r^, 

b.  Depth  z = 0.5  r^, 

c.  Depth  z = 0.9  r^ 

Fig.  16.  Absorbed  dose  along  the  central  axis  of  a square  fields  with  sides  of  1,  2,  4,  and  8 mm, 
for  a 160-MeV  beam.  The  dotted  curve  is  for  a laterally  unbounded  broad  beam. 

Fig.  17.  Absorbed  dose  along  the  central  axis  of  rectangular  fields  that  are  imbounded  in  one 

dimension  and  have  widths  of  1,  2,  4,  and  8 mm  in  the  other  dimension,  for  a 160-MeV 
beam.  The  dotted  curve  is  for  a laterally  unbounded  broad  beam. 

Fig.  18.  Reduction  factors  and  absorbed-dose  values  for  a rectangular  1 mm  x 10  mm  field,  for  a 
160-MeV  beam. 

a.  Central-axis  depth  dose 

b.  At  depth  z = 0.5  r^,,  as  a function  of  radial  distance  along  a line  that  starts  at 
the  center  of  the  field  and  is  parallel  to  the  long  side  of  the  field 

c.  Similar  to  b,  but  for  a line  that  is  parallel  to  the  short  side  of  the  field. 
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Fig. 5 
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Fig. 6 
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Fig. 7a 
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Fig. 7b 
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Fig.  8 
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Fig. 9b 
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Fig. 14c 
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Fig. 15a 
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Fig. 15c 
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Fig. 16 
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Fig. 17 
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Fig. 18a 
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Fig. 18b 
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Fig. 18c 


Appendix  1.  Files  on  Disk 

Fortran  source  codes,  executable  codes  (for  use  on  IBM-compatible  personal  computers), 
data  files  and  sample  output  files  are  stored  on  a 3.5  inch  1.44  Mb  floppy  disk.  Included  are  the 
following  files: 


AXPLOT.EXE 

FREC.EXE 

MOLC2.COF 

PTGPOL.FOR 

AXPLOT.FOR 

FREC.FOR 

MORAD.160 

PTPOL.EXE 

COMPOS.WAT 

FRECl.l 

MORAD.EXE 

PTPOL.FOR 

FCIR.l 

FREC1.EXE 

MORAD.FOR 

PTRAD.160 

FCIR.2 

FRECl.FOR 

MORAD10.160 

PTRAD.EXE 

FCIR.EXE 

FREC2.EXE 

MORAD160.050 

PTRAD.FOR 

FCIR.FOR 

FREC2.FOR 

MORAD50.160 

RADPLOT.EXE 

FREC.l 

MCORR.OOl 

MORAD90.16 

RADPLOT.FOR 

FREC.2A 

MCORR.008 

PDEPTH.ARR 

STOPRANG.WAT 

FREC.2B 

MOLCl.COF 

PTGPOL.EXE 

ZRFILE 
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Appendix  2.  Source  Code  Listings 


The  following  programs  are  listed: 


PTPOL 

Energy  deposition  distributions  from  monoenergetic  beams 

PTGPOL 

Energy  deposition  distributions  from  beams  with  a Gaussian  spectrum 

MORAD 

Radial  dose  distributions  from  Moli^re’s  theory 

PTRAD 

Radial  dose  distributions  from  Monte  Carlo  progrmi  PTRAN 

FCIR 

Calculation  of  3-D  dose  distributions  from  circular  fields 

FREC 

Calculation  of  3-D  dose  distribution  from  rectangular  field,  along  a 
line  starting  at  the  center  of  the  field. 

FRECl 

Calculation  of  3-D  dose  distribution  from  rectangular  field,  for  a set 
of  (x,y)  positions 

FREC2 

Same  as  FRECl,  using  double  numerical  quadrature 

AXPLOT 

Plotting  of  central-axis  absorbed-dose  distributions 

RADPLOT 

Plotting  of  radial  absorbed-dose  distributions 
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PROGRAM  PTPOL 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


10  May  93.  Written  by  Martin  J.  Berger,  NIST. 

Calculates  energy  deposition  per  unit  depth,  as 
function  of  depth,  for  a monoenergetic  beam  with 
with  an  energy  between  SO  and  250  MeV.  Results  are 
obtained  by  Interpolation  In  a database  of  processed 
results  obtained  with  the  Monte  Carlo  program  PTRANIO 
at  25  energies  between  250  MeV  and  50  MeV. 

Required  Input  file: 

PDffTH.ARR,  processed  depth-dose  distributions 

The  following  subroutines  are  called: 

SCOFD  Calculates  cubic  spline  coefficients 

BSPOLD  Evaluates  cubic  spline 

BSP0L02  Evaluetes  first  and  second  derivative 

of  cubic  spline 

SPMAX  Finds  maximum  of  cubic-spline  function 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-2) 

DIMENSION  T(30),TL(30),Z(50),D(30,50),DL(30), 

1 AX1(30),BX1(30),CX1(30),DX1(30), 

5 AR(30),BR(30),CR(30),OR(30), 

6 RG(3O),RGL(3O),OINT(5O),OINTl(50) 

CHARACTER  0UTPUT*30,LINE*80 

DATA  EPS/1. OD-05/ 

1 FORMAT (IH  ) 

5 FORMAT(A) 

PRINT  Enter  energy  (MeV):  ' 

READ  *,  EBEAM 

IF(EBEAM.LE.250.0DO.AND.EBEAM.GE.50.0DO)  GO  TO  7 
PRINT  6 

6 F0RMAT('  Energy  Is  out  of  range.') 

STOP 

7 PRINT  Enter  name  of  output  file:  ' 

READ  5,  OUTPUT 

OPEN  (8, OUTPUT) 

WRITE  (8,8)  OUTPUT 

8 FORMAT ( 'Program  PTPOL,  output  file  ',A) 

WRITE  (8,1) 

OPEN  (7,'PDEPTH.ARR') 

READ  (7,5)  LINE 
READ  (7,«)  NNAX,LMAX 
READ  (7,*)  (Z(N),N-1,NMAX) 

DO  SO  L-1,LMAX 
READ  (7,*)  T(L),RG(L) 

TL(L)-LOG(T(L)) 

RGL(L)-LOG(RG(L)) 

READ  (7,*)  (D(L,N),N-1,NMAX) 

DO  40  N-1,NHAX 

IF(D(L,N).LE.O.ODO)  0(L,N)-1 . OD-08 
40  CONTINUE 
50  CONTINUE 

CALL  SCOFD(TL,RGL,LMAX,AR,BR,CR,DR) 

EBEAML-LOG( EBEAM) 

CALL  BSPOLD(EBEAML,TL,AR,BR,CR,DR,LMAX,RES) 

RANGE-EXP (RES) 

DO  80  N-1,NMAX 
DO  70  L-1,LMAX 
0L(L)-L0G(D(L,N)) 

70  CONTINUE 

CALL  SCOFD(TL,OL,LHAX,AX1,BX1,CX1,OX1) 

CALL  BSPOLD ( EBEAML , TL , AXl , BXl , CXI , DXl , LMAX , RES ) 

DINT(N)-EXP(RES) 

80  CONTINUE 

CALL  SPMAX(Z,DINT,NNAX,EPS,ZH,DMAX) 

DO  100  N-1,NHAX 
100  0INT1(N)-DINT(N)/DHAX 

110  WRITE  (8,111) 

111  FORMAT('  EBEAM  - average  beam  energy,  MeV') 

WRITE  (8,112) 

112  FORMAT('  RANGE  - CSDA  range  at  energy  EBEAM,  g/cra2') 

WRITE  (8,113) 

113  F0RHAT('  ZM  - depth  (In  units  of  RANGE)  at  which  energy-deposi 
Itlon  curve  peaks' ) 

WRITE  (8,114) 

114  FORNAT('  DHAX  - energy  deposition  per  unit  depth  at  peak,  MeV  cm 

12/g') 

WRITE  (8,115) 

115  F0RHAT('  NMAX  - number  of  depths') 

WRITE  (8,1) 

WRITE  (8,116) 

116  F0RMAT(4X, ' EBEAM' ,4X, 'RANGE' ,7X, 'ZH' ,5X, 'DHAX' ,3X,  'NMAX' ) 

WRITE  (8,117)  EBEAM, RANGE, ZM, DHAX, NMAX 

117  FORMAT(F9.2,3F9.4,I7) 

WRITE  (8,1) 

WRITE  (8,120) 

120  FORHAT( 'Depths,  In  units  of  RANGE') 

WRITE  (8,130)  (Z(N),N-1,NMAX) 

130  FORMAT(6F12.3) 

WRITE  (8,140) 

140  FORMAT('dO/dz,  energy -deposit Ion  distribution,  MeV  Cffl2/g,  per  Inci 
Ident  proton' ) 


WRITE  (8,150)  (DINT(N),N-1,NMAX) 

150  FORMAT (6F 12. 6) 

WRITE  (8,160) 

160  FORHAT( 'Relative  energy -deposit Ion  distribution  (peak  value  unity) 
1') 

WRITE  (8,150)  (0INT1(N),N-1,NMAX) 

STOP 

END 

SUBROUTINE  SPMAX(X,Y,NMAX,EPS,XM,YMAX) 

Z 3 Mar  93.  Finds  maximum  of  spline  function. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  X(1000),Y(1000),A(1000),B(1000),C(1000),D(1000) 

YMAX-0.0 

XM-X(l) 

DO  20  N-1,NMAX 
IF(Y(N)-YMAX)20,20,10 
10  YMAX-Y(N) 

XM-X(N) 

20  CONTINUE 

CALL  SC0FD(X,Y,NMAX,A,B,C,0) 

30  CALL  BSPOL02(XM,X,A,B,C,0,NHAX,RES,RES1,RES2) 

XNNEW-XM-RES1/RES2 
RAT-ABS ( ( XMNEW-XM ) /XM ) 

IF(RAT.LT.EPS)  GO  TO  40 
XM-XMNEW 
GO  TO  30 
40  YMAX-RES 
RETURN 
END 

SUBROUTINE  BSPOLD2(S,X,A,B,C,0,N,G,G1,G2) 

C 2 Mar  93.  Evaluates  first  and  second  derivative  of  cubic-spline 
C function. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  X(IOOO) ,A(1000) ,B(1000) ,C(1000) ,0(1000) 

IF  (X(l).GT.X(N))  GO  TO  10 
IDIR-0 
MLB-0 
MUB-N 
GO  TO  20 
10  IDIR-1 
MLB-N 
MUB-O 

20  IF  (S.GE.X(MUB+IDIR))  GO  TO  60 
IF  (S.LE.X(NLB+1-IDIR))  GO  TO  70 
ML-MLB 
MU-MUB 
GO  TO  40 

30  IF  (lABS(MU-NL).LE.l)  GO  TO  80 
40  HAV-(ML^HU)/2 

IF  (S.LT.X(MAV))  GO  TO  50 
ML-MAV 
GO  TO  30 
50  MU-MAV 
GO  TO  30 

60  MU-MUB4^2*IDIR-1 
GO  TO  90 

70  MU-MLB-2*IDIR+1 
GO  TO  90 
80  MU-MU^IDIR-l 
90  Q-S-X(NU) 

G-( (D(MU)*Q+C(MU) )*Q+B(MU) )*Q+A(MU) 

G1-B(MU)+2.0*C(MU)*Q+3.0*D(MU)*Q«Q 

G2-2.0*C(MU)+6.0*D(MU)*Q 

RETURN 

EW) 

SUBROUTINE  SC0FD(X, F,NHAX,A,B,C,0) 

C 1 JAN  87.  DOUBLE  PRECISION 
C IF  S LIES  BETWEEN  X(N)  AND  X(M^l),  THEN 
C F(S)-(  (0(M)*S<-C(M)  )«S+B(M)  )«S+A(M) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  X(IOOO) ,F(1000) ,A(1000) ,B(1000) ,C(1000) ,D(1000) 

Ml-2 

M2-W<AX-1 

S-0.0 

DO  10  N-1,H2 
D(M)-X(M+1)-X(M) 

R-(F(M+1)-F(M))/D(H) 

C(H)-R-S 
10  S-R 
S-0.0 
R-0.0 
C(l)-0.0 
C(NNAX)-0.0 
DO  20  H-H1,H2 
C(H)-C(M)+R*C(M-1) 

B(M)-(X(M-1)-X(M+1))*2.0D0-R*S 

S-D(H) 

20  R-S/B(N) 

NR-N2 

DO  30  N-M1,M2 

C(MR)-(0(MR)*C(MR+1)-C(HR))/B(MR) 

30  HR-MR-1 

DO  40  M-1,M2 
S-D(N) 

R-C(M+1)-C(M) 

D(N)-R/S 
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C(H)-C(M)*3.000 

B(M)-(F(M+1)-F(H))/S-(C(H)+R)*S 
40  A(M)-F(M) 

RETURN 

END 

SUBROUTINE  BSPOLD(S,X,A,B,C,0,N,G) 

C 1 JAN  87.  DOUBLE  PRECISION 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  X(IOOO) ,A(1000) ,B(1000) ,C(1000) ,0(1000) 
IF  (X(l).GT.X(N))  GO  TO  10 
IDIR-0 
MLB-O 
HUB-N 
GO  TO  20 
10  IDIR-l 
MLB-N 
MUB-O 

20  IF  (S.GE.X(MUB^IOIR))  GO  TO  60 
IF  (S.LE.X(MLB+1-IDIR))  GO  TO  70 
HL-NLB 
HU-NUB 
GO  TO  40 

30  IF  (lABS(MU-NL).LE.l)  GO  TO  80 
40  HAV-(HL+HU)/2 

IF  (S.LT.X(MAV))  GO  TO  50 
ML-MAV 
GO  TO  30 
50  HU-MAV 
GO  TO  30 

60  NU-MUB+2« IDIR-l 
GO  TO  90 

70  NU-NLB-2*I0IR^1 
GO  TO  90 
80  MU-MU^IDIR-l 
90  Q-S-X(HU) 

6.((0(HU)*Q+C(MU))*Q+B(MU))*Q+A(MU) 

RETURN 

END 


) 
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PROGRAM  PTGPOL 

10  May  93.  Written  by  Martin  J.  Berger,  NIST. 

Calculates  energy  deposition  per  unit  depth,  as 
function  of  depth,  for  bean  with  a Gaussian  spectrum. 
The  spectrum  must  be  specified  In  terms  of  a mean 
energy,  and  a relative  standard  deviation  (expressed 
as  a percentage  of  the  mean  energy.  PTGPOL  relies 
on  Interpolation  In  a database  of  processed  results 
obtained  with  Monte  Carlo  program  PTRANID  for  proton 
beans  with  25  energies  between  250  MeV  and  50  MeV. 

Required  input  file: 

POEPTH.ARR,  processed  depth-dose  distributions 

The  following  subroutines  are  called; 

SCOFO  Calculates  cubic  spline  coefficients 

BSPOLO  Evaluates  cubic  spline 

BSP0L02  Evaluetes  first  and  second  derivative 

of  cubic  spline 

SPHAX  Finds  maximum  of  cubic-spline  function 

IMPLICIT  DOUBLE  PRECISION  <A-H,0-2) 

DIMENSION  T(30), TL(30),Z(50),D(30, 50), DL(30), 00(201), DOREL(201), 

1 E(401),WT(401),RGINT(401),AX1(30),BX1(30),CX1(30),DX1(30), 

2 AX3(50),BX3(SO),CX3(50), 0X3(50), AR(30),BR(30),CR(30),0R(30), 

3 RG(30),RGL(30),OIN(50),DINT(401,50),OS(401,52),2REF(201) 

CHARACTER  INPUT*16,OUTPUT*16,LINE*80 

DATA  EPS/1. OD-05/ 

I FORMAT (IH  ) 

5 FORMAT(A) 

PRINT  Enter  average  beam  energy  (MeV):  ' 

READ  *,  EAV 

PRINT  Enter  relative  standard  deviation  (percent):  ' 

READ  *,  P 

SIGMA-0. 01*P*EAV 

PRINT  Suggested  values:  4.0' 

PRINT  Enter  cut-off  for  Gaussian  (In  units  of  standard  devlatl 
Ion:  ' 

READ  •,  CUT 
ENMAX-EAV+CUT*SI6HA 
IF(ENHAX.LE.2S0.0D0)  GO  TO  6 

STOP  ' STOP:  Some  energies  In  spectrum  are  above  250  MeV.' 

6 ENHIN-EAV-CUT*SIGMA 
IF(ENHIN.GE.50.0D0)  GO  TO  7 

STOP  ' STOP:  Some  energies  In  spectrum  are  below  50  MeV.' 

7 PRINT  *,'  Suggested  value:  401' 

PRINT  *,'  Enter  number  of  energies  In  Gaussian  that  are  sampled:  ' 
READ  *,  KMAX 

PRINT  *,'  Suggested  choice:  ZRFILE  (default)' 

PRINT  *,'  Enter  name  of  file  with  list  of  depths:  ' 

READ  5,  INPUT 

PRINT  Enter  name  of  output  file:  ' 

READ  5,  OUTPUT 
OPEN  (8, OUTPUT) 

WRITE  (8,10)  OUTPUT 

10  FORMAT( ' Program  PTGPOL,  output  file  ',A) 

WRITE  (8,11)  INPUT 

II  FORMAT  ('List  of  depths  from  file  ',A) 

WRITE  (8,1) 

OPEN  (7, INPUT) 

READ  (7,«)  NRMAX 

READ  (7,*)  (ZREF(N),N-1,NRHAX) 

CLOSE  (7) 

EHIN-EAV-CUT*SIGMA 
EHAX-EAV-fCUT*SIGNA 
OE- ( EMAX- EM IN ) /REAL (KMAX- 1 ) 

00  15  K-1,KHAX 
E(K)-EMIN40E*REAL(K-1) 

WT{K)-0.0 

IF(E(K).GT.EHAX)  GO  TO  15 
IF(E(K).LT.EMIN)  GO  TO  15 
ARG-0.5*((E(K)-EAV)/SI6MA)**2 
WT(K)-EXP(-ARG) 

15  CONTINUE 
SlM-0.0 
DO  20  K-1,KHAX 
20  SUH-SUM^WT(K) 

DO  30  K-1,KMAX 
30  WT(K)-WT(K)/SUM 

OPEN  (7, 'POEPTH.ARR') 

READ  (7,5)  LINE 
READ  (7,«)  NHAX,LHAX 
READ  (7,*)  (Z(N),N-1,NMAX) 

DO  50  L-1,LMAX 
READ  (7,*)  T(L),RG(L) 

TL{L)-LOG(T(L)) 

RGL(L)-LOG(RG(L)) 

READ  (7,*)  (D(L,N),N-1,NHAX) 

DO  40  N-1,NHAX 

IF(D(L,N).LE. 0.000)  D(L,N)-1. 00-08 
40  CONTINUE 
50  CONTINUE 

CALL  SCOFO(TL,RGL,LNAX,AR,BR,CR,OR) 


EAVL-LOG(EAV) 

CALL  BSPOLD(EAVL,TL,AR,BR,CR,OR,LMAX,RES) 

RGREF-EXP(RES) 

DO  60  K-1,KMAX 
EL-LOG(E(K)) 

CALL  BSPOLO(EL,TL,AR,BR,CR,OR,LMAX,RES) 

60  RGINT(K)-EXP(RES) 

DO  90  N-1,NMAX 
DO  70  L-1,LHAX 
DL(L)-L0G(0(L,N)) 

70  CONTINUE 

CALL  SCOFD( TL , DL , LMAX , AXl , BXl , CXI , DXl ) 

DO  80  K-1,KNAX 
EL-LOG(E(K)) 

CALL  BSPOLO ( EL , TL , AXl , BXl , CXI , DXl , LHAX, RES ) 

DINT(K,N)-RES 
80  CONTINUE 
90  CONTINUE 

DO  140  K-1,KMAX 
DO  100  N-1,NNAX 
DIN(N)-DINT(K,N) 

100  CONTINUE 

CALL  SCOFD( Z , DIN, NMAX , AX3 , BX3 , CX3 , DX3 ) 

DO  130  N-1, NRMAX 
0D(N)-0.0 

Z INT-ZREF ( N ) ♦ ( RGREF/RG I NT ( K ) ) 

IF(ZINT-1. 04)120, 110, 110 
110  0S(K,N)-0.0 
GO  TO  130 

120  CALL  BSPOLO(ZINT,Z,AX3,BX3,CX3, 0X3, NMAX, RESl) 

0S(K,N)-EXP(RES1) 

130  CONTINUE 
140  CONTINUE 

DO  160  N-1, NRMAX 
DD(N)-0.0 
DO  150  K-1,KMAX 
DD(N)-DD(N)+WT(K)*DS(K,N) 

150  CONTINUE 
160  CONTINUE 

CALL  SPNAX(ZREF,OD, NRMAX, EPS, ZM.OHAX) 

DO  180  N-1, NRMAX 
180  OOREL(N)-DO(N)/DHAX 
WRITE  (8,191) 

191  FORMAT('  EAV  - average  beam  energy,  MeV') 

WRITE  (8,192) 

192  F0RMAT('  SIGMA  - standard  deviation  of  beam  spectrum,  MeV') 

WRITE  (8,193) 

193  FORMAT('  P - 100*SIGMA/EAV,  relative  standard  deviation') 

WRITE  (8,194) 

194  FORHAT('  CUT  - cut-off  for  Gaussian  beam  (In  units  of  SIGMA)') 
WRITE  (8,195) 

195  FORHAT('  RGREF  - CSDA  range  at  energy  EAV') 

WRITE  (8,196) 

196  FORHAT('  ZM  - depth  (In  units  of  RANGE)  at  which  energy-depos 
lltlon  curve  peaks') 

WRITE  (8,197) 

197  FORHAT('  DHAX  - energy  deposition  per  unit  depth  at  peak,  MeV  cm 
12/g') 

WRITE  (8,198) 

198  FORMAT('  NRMAX  - number  of  depths') 

WRITE  (8,1) 

WRITE  (8,220) 

220  FORMAT (6X, 'EAV' ,4X, 'SIGMA' ,8X, 'P' ,6X, 'CUT' ,4X, 'RGREF' ,7X, 'ZM' , 

1 5X, 'DHAX', 3X, 'NRMAX') 

WRITE  (8,230)  EAV, SIGMA, P, CUT, RGREF, ZM, DHAX, NRMAX 
230  F0RHAT(4F9.2,3F9.4,I8) 

WRITE  (8,1) 

WRITE  (8,240) 

240  FORNAT( 'Depths,  In  units  of  RGREF') 

WRITE  (8,250)  (ZREF(N) , N-1, NRMAX) 

250  FORMAT (6F12. 6) 

WRITE  (8,260) 

260  FORNAT('dD/dz,  energy-deposition  distribution,  MeV  cmZ/g,  per  Inci 
Ident  proton' ) 

WRITE  (8,250)  (00(N) , N-1, NRMAX) 

WRITE  (8,270) 

270  F0RMAT( 'Relative  energy-deposition  distribution  (peak  value  unity' 
1) 

WRITE  (8,250)  (DOREL(N) , N-1 , NRMAX) 

STOP 

END 

SUBROUTINE  SPMAX(X,Y,NHAX,EPS,XM,YMAX) 

3 Mar  93.  Finds  ouxlmin  of  cubic-spline  function. 

IW>LICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  X(IOOO) ,Y(1000) ,A(1000) ,B(1000) ,C(1000) ,0(1000) 

YNAX-0.0 

XN-X(l) 

DO  20  N-1, NMAX 
IF(Y(N)-YMAX)20,20,10 
10  YMAX-Y(N) 

XM-X(N) 

20  CONTINUE 

CALL  SCOFO(X,Y,NMAX,A,B,C,D) 

30  CALL  BSP0L02(XM,X,A,B,C,0,NHAX,RES,RES1,RES2) 

XNNEW-XN-RES1/RES2 
RAT-ABS ( ( XHNEW-XM ) /XM ) 
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IF(RAT.LT.EPS)  GO  TO  40 
XH-XMNEW 
GO  TO  30 
40  YNAX-RES 
RETURN 
END 

SUBROUTINE  BSPOL02(S,X,A,B,C,D,N,G,G1,G2) 

C 2 Har  93.  Evaluates  Hrst  and  second  derivative. 

C of  cubic-spline  function 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  X(IOOO) ,A(1000) ,B(1000) ,C(1000) ,D(1000) 

IF  (X(l).GT.X(N))  GO  TO  10 
IDIR-0 
MLB-O 
MUB-N 
GO  TO  20 
10  IDIR-1 
MLB-N 
MUB-0 

20  IF  (S.GE.X(HUB^IDIR))  GO  TO  60 
IF  (S.LE.X(NLB4^1-I0IR)}  go  to  70 
ML-HLB 
MU-MUB 
GO  TO  40 

30  IF  (lABS(HU-ML).LE.l)  GO  TO  80 
4C  Kr.v-(ML+MU)/2 

!S.LT.X(HAV))  go  to  50 

to  30 
50  NU-HAV 
GO  TO  30 

60  HU-HUB4^2*I0IR-1 
GO  TO  90 

70  HU-HLB-2*IDIR+1 
GO  TO  90 
80  NU-NU4^I0IR-1 
90  Q-S-X(MU) 

6-( (D(MU)*Q+C(MU) )*Q+B(MU) )*Q+A(HU) 

Gl-B(MU)+2. 0*C(MU)*Q+3. 0*D(HU)*Q*Q 

G2-2.0«C(MU)4-6.0«D(HU)*q 

RETURN 

END 

SUBROUTINE  SCOFO(X,F,NNAX,A,B,C,D) 

C 1 JAN  87.  DOUBLE  PRECISION 
C IF  S LIES  BETWEEN  X(N)  AND  X(N+1),  THEN 
C F(S).(  (D(M)*S+C(H)  )*S4^B(N)  )*S+A(M) 

II^LICIT  DOlffiLE  PRECISION  (A-H,0-Z) 

DIMENSION  X(1000),F(1000),A(1000),B(1000),C(1000),D(1000) 
Nl-2 

N2-NHAX-1 

S-0.0 

DO  10  H-1,H2 
D(N)-X(M>1)-X(N) 

R-(F(M+1)-F(M))/D(M) 

C(H)-R-S 
10  S-R 
S-0.0 
R-0.0 
C(l)-0.0 
C(NMAX)-0.0 
DO  20  M-M1,H2 
C{H)-C(H)+R*C(M-1) 

B(H)-(X(M-1)-X(M+1))*2.0D0-R*S 

S-O(N) 

20  R-S/B(N) 

NR-N2 

DO  30  N-M1,H2 

C(MR)-(D(MR)*C(MR+l)-C(MR))/B(HR) 

30  NR-HR-1 
DO  40  N-l,H2 
S-O(N) 

R-C(M+l)-C(M) 

0(N)-R/S 

C(N)-C(N)«3.0D0 

B{M)-(F(M+1)-F(M))/S-(C(N)+R)*S 
40  A(N)-F(H) 

RETURN 

END 

SUBROUTINE  BSPOLO(S,X,A,B,C,D,N,G) 

C 1 JAN  87.  DOUBLE  PRECISION 

IMPLICIT  DOUBLE  PRECISION  <A-H,0-Z) 

DIMENSION  X(IOOO) ,A(1000) ,B(1000) ,C(1000) ,D(1000) 

IF  (X(1).6T.X(N))  GO  TO  10 
IDIR-O 
MLB-O 
MUB-N 
GO  TO  20 
10  IDIR-1 
MLB-N 
MUB-0 

20  IF  (S.GE.X(MUB+IDIR))  GO  TO  60 
IF  (S.LE.X(NLB4^1-IDIR))  go  to  70 
ML-MLB 
MU-MUB 
GO  TO  40 

30  IF  (lABS(MU-ML).LE.l)  GO  TO  80 


40  HAV-(NL+MU)/2 

IF  (S.LT.X(MAV))  GO  TO  50 
ML-MAV 
GO  TO  30 
50  MU-MAV 
GO  TO  30 

60  MU-MUB+2*I0IR-1 
GO  TO  90 

70  MU-MLB-2«IDIR4^1 
GO  TO  90 
80  MU-MU-^IDIR-1 
90  Q-S-X(MU) 

G-( (D(MU)*Q+C(MU) )*Q+B(MU) )*Q+A(MU) 

RETURN 

END 
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PROGRAM  HORAO 
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8 July  93.  Written  by  Martin  J.  Berger,  NIST. 

Evaluates  radial  multiple  scattering  distribution 
f(rho)  of  Mollere  (Z.Naturf. 10a, 177-211, 1955),  as 
function  of  the  radial  distance  rho. 


The  expansion  coefficients  In  the  Mollere  formula, 
are  supplied  by  files  MOLCl.COF  and  M0LC2.C0F. 


For  application  to  a water  medium,  the  following 
data  files  are  required  as  Input: 


STOPRANG.WAT  Proton  stopping  powers  and  ranges 
In  water 


COMPOS. WAT  Composition  data  for  water 

MCORR.OOl  and  Correction  factors  for  H and  0 that 
HCORR.008  convert  Thomas-Fermi  to  Hartree- 

Fock  screening 

The  following  subroutines  are  used: 


SCOFD:  calculates  cubic  spline  coefficients 

BSPOLO:  evaluates  cubic  spline 

GRALO:  carries  out  numerical  quadrature 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-2) 

COMMON  X(401) ,0(401 ) , COF1(401 ) ,COF2(401 ) , 

1 DOR(3001),MZ(20),WT(20),CP(20), 

2 APL(20),ALFB(20),FAC1(20),FAC8(20), 

3 AM1(20),BH1(20),CH1(20),DM1{20), 

4 AN8(20),BM8(20),CM8(20),DM8(20) 

COMMON  THAI ( 451 ) , THBl ( 451 ) , THA2 ( 451 ) , THB2 ( 451 ) , 

1 FA1{451),FB1(451),FA2(451),FB2(451), 

2 AS1(451),BS1(451),CS1(451),0S1(451), 

3 AL1(451),BL1(451),CL1(451),DL1{451), 

4 AS2(451),BS2(451),CS2(451),DS2(451), 

5 AL2(451),BL2(451),CL2(451),DL2(451) 

COMMON  ENG( 133) , ENGL( 133) ,STP (133 ) , RANG ( 133 ) , RANGL( 133) , 

1 A1 (574) , B1 (574) , Cl (574) , 01(574), 

2 A2(574),B2(574),C2(574), 02(574), 

3 A3(101),B3(101),C3(101), 03(101), 

3 S(1001),SFAC(1001),TAUFAC(1001),TAU2(1001),BETQ(1001),BET(1001), 

4 aVVN0(  1001  ),ZR(  120) 

DIMENSION  UIN(IOO) 

CHARACTER*72  HEAD1,HEAD2,TA6 
CHARACTER*16  INPUT, OUTPUT 

DATA  UIN/3. 6, 4. 0,4. 3, 4. 6, 4. 8, 5. 0,5*5. 1,3*5. 2, 3*5. 3, 3*5. 4, 

1 4*5.5,5*5.6,5*5.7,6*5.8,7*5.9,8*6.0,10*6.1,11*6.2,12*6.3, 

2 12*6.4/ 

DATA  FCON/1130.0DO/,LHAX/1001/ 

DATA  AVOG/6 . 02213670+23/ , EMASS/O . 51 09990600/ , PHASS/938. 2723100/ , 

1 FINE/137 . 03598900/ , RCLASS/2 . 817940920- 13/ , EUL/0 . 577215664900/ 
DATA  BEPS/1 . OE-09/, ZRMAX/0 . 98/, THHAX/10 . 000/ , NMAX/99/ , IFAN/1/ 

1 FORMAT (IH  ) 

5 FORHAT(A) 

PI-4. 000* AT AN (1.000) 

RATMQ-(EMASS/PMASS)**2 

CON-4 . ODO*PI*AVOG* ( RCLASS**2 ) *RATMq 

TFC0N-( (3. ODO*PI)**(2. 000/3. 000) )/(2.000**(7. 000/3. 000)) 

C0N1-RATMQ/(TFC0N*FINE)**2 

PRINT  *,'  Enter  Initial  proton  energy  (MeV):  ' 

READ  *,  TIN 

PRINT  *,'  Path  lengths  are  to  be  expressed  in  units  of  CSDA  range' 
PRINT  *,'  Options  for  path  length  Input:' 

PRINT  *,'  1)  From  file' 

PRINT  *,'  2)  Specify  set  of  path  lengths  In  terms  of  maximum' 

PRINT  *,'  value,  ZRMAX,  and  number  of  values,  NMAX.' 

PRINT  *,'  3)  Use  default  values  ZRMAX-0.98,  NHAX-98.' 

PRINT  *,'  Enter  choice:  ' 

READ  *,  INF 
GO  TO  (6,7,8),  INF 

6 PRINT  *,'  Enter  name  of  path- length  Input  file:  ' 

READ  5,  INPUT 

OPEN  (7, INPUT) 

READ  (7,*)  NMAX 

READ  (7,*)  (ZR(N),N-1,NMAX) 

CLOSE  (7) 

GO  TO  12 

7 PRINT  *,'  Enter  maximum  path  length,  ZRMAX:  ' 

READ  «,  ZRMAX 

PRINT  *,'  Enter  nuaber  of  path  lengths  (no  greater  than  120):  ' 
READ  *,  NMAX 
GO  TO  9 

8 ZRHAX-0.98 
NMAX-98 

9 DZR-ZRNAX/OBLE(NMAX) 

00  10  N-1,NHAX 

10  ZR(N)-OZR*DBLE(N) 

12  PRINT  *,'  Enter  nianber  of  radial  distances  (no  greater  than  401): 
1' 

READ  *,  KMAX 

PRINT  *,'  Suggested  name:  MORAO.nnn' 


PRINT  *,'  Enter  name  of  output  file:  ' 

READ  5,  OUTPUT 
OPEN  (8, OUTPUT) 

WRITE  (8,13)  OUTPUT 

13  FORMAT( 'Program  MORAD,  output  file  ',A) 
WRITE  (8,1) 

IF(NMAX.EQ.l)  WRITE  (8,14)  NMAX 

14  F0RMAT(I6,'  depth') 

IF(WfAX.NE.l)  WRITE  (8,15)  NMAX 

15  F0RMAT(I6,'  depths') 

IF(KMAX.EQ.l)  WRITE  (8,16)  KMAX 

16  FORMAT(I6,'  radial  distance') 

IF(KMAX.NE.l)  WRITE  (8,17)  KMAX 

17  F0RMAT(I6,'  radial  distances') 

WRITE  (8,1) 

OPEN  (UNIT-7, FILE-' STOPRANG.WAT') 

READ  (7,*)  JMAX 

READ  (7,*)  ( ENG (J),J-1, JMAX) 

READ  (7,*)  (STP(J),J-1,JHAX) 

READ  (7,*)  (RANG(J),J-1,JHAX) 

CLOSE  (7) 

DO  18  J-1,JNAX 
ENGL(J)-LOG(ENG(J)) 

18  RANGL(J)-LOG(RANG(J)+l. 00-09) 

CALL  SCOFD(ENGL,RANGL,JMAX,A1,B1,C1,01) 

CALL  SC0FD(RANGL,ENGL,JMAX,A2,B2,C2,02) 

OPEN  (7, 'COMPOS. WAT') 

READ  (7,5)  MAT 
READ  (7,*)  MMAX 

READ  (7,20)  (M2(M),WT(H),M-1,MMAX) 

20  F0RMAT(6(I3,F9.6)) 

CLOSE  (7) 

TINL-LOG(TIN) 

CALL  BSPOLO(TINL, ENGL,A1, B1 ,C1, 01 , JMAX, RES) 
R6IN-EXP(RES) 

OPEN  (UNIT-7, FILE-'HOLCl.COF') 

READ  (7,5)  HEADl 

READ  (7,5)  TAG 

READ  (7,5)  TAG 

READ  (7,5)  TAG 

READ  (7,*)  JMA1,JMB1,THCUT1 

READ  (7,*)  (THA1(J),J-1,JMA1) 

READ  (7,*)  (FA1(J),J-1,JMA1) 

READ  (7,*)  (THB1(J),J-1,JMB1) 

READ  (7,*)  (FB1(J),J-1,JMB1) 

CLOSE  (7) 

CALL  SCOFO(THA1,FA1,JMA1,AS1,BS1,CS1,OS1) 
CALL  SC0F0(THB1 , FBI , JMBl , ALl , BLl , CLl , DLl ) 
OPEN  (UNIT-7, FILE-' M0LC2.C0F') 

READ  (7,5)  HEA02 

READ  (7,5)  TAG 

READ  (7,5)  TAG 

READ  (7,5)  TAG 

READ  (7,*)  JMA2,JMB2,THCUT2 

READ  (7,*)  (THA2(J),J-1,JMA2) 

READ  (7,*)  (FA2(J),J-1,JMA2) 

READ  (7,*)  (THB2(J),J-1,JMB2) 

READ  (7,*)  (FB2(J),J-1,JHB2) 

CLOSE  (7) 

CALL  SCOFD(THA2,FA2,JMA2,AS2,BS2,CS2,DS2) 
CALL  SC0FD( THB2 , FB2 , JNB2 , AL2 , BL2, CL2 , DL2 ) 
OPEN  (7, 'MCORR.OOl' ) 

READ  (7,5)  TAG 

READ  (7,*)  MCMAX 

DO  30  MC-l, MCMAX 

READ  (7,*)  ALFB(HC),FAC1(MC) 

30  CONTINUE 
CLOSE  (7) 

CALL  SC0FD( ALFB , FACl , MCMAX, AMI , BMl , CMl , DM1 ) 
OPEN  (7,'MC0RR.008') 

READ  (7,5)  TAG 

READ  (7,*)  MCMAX 

DO  40  MC-l, MCMAX 

READ  (7,*)  ALFB(MC),FAC8(HC) 

40  CONTINUE 
CLOSE  (7) 

CALL  SCOFD( ALFB, FAC8, MCMAX, AM8,BH8,CN8, DM8} 

DO  280  N-1,NMAX 

RFRAC-ZR(N) 

PATH-RGIN«RFRAC 

0S-PATH/REAL(LMAX-1) 

DELTA-OS/3.000 

DO  50  L-1,LHAX 

S(L)-0S*REAL(L-1) 

SFAC(L)-(1.000-S(L)/PATH)**2 

R6L-LOG(RGIN-S(L)) 

CALL  BSP0L0(RGL,RANGL,A2,B2,C2,D2,JHAX,RES) 
E-EXP(RES) 

TAU-E/PMASS 

TAUFAC( L) -( (TAU+1 . 0 ) /TAU/ (TAU+2 . 0 ) ) **2 
TAU2{L)-TAU*(TAU+2.0D0) 
BETQ(L)-TAU*(TAU+2.0)/( (TAU+1. 0)**2) 

50  BET(L)-SQRT(BETg(L)) 

DO  80  H-1,MHAX 
IZ-MZ(H) 

Z-REAL(IZ) 
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IF(IZ.EQ.l)  ATVT-1.0079D0 
IF(IZ.EQ.8)  ATWT-15. 999600 
ZTT-Z**( 2. 000/3.000) 

ZFT-ZTT**2 

UINF-UIN(IZ) 

CONZ-CON* ( Z*Z/ATWT ) *WT{ H ) 

00  60  L-1,LHAX 

60  GRANO(L)-TAUFAC(L)*SFAC(L) 

CALL  GRAL0(0ELTA,GRAN0,LHAX,RES1) 

CP(H)-C0NZ*RES1 
DO  70  L-1,LHAX 
ALPH-Z/FINE/BET(L) 

FAC-1 . 1300+3. 76D0*(ALPH**2) 

IF(H.EQ.l)  CALL  BSP0LD(ALPH,ALFB,AH1,BH1,CH1,DH1,HCHAX,ANS) 
IF(H.EQ.Z)  CALL  BSP0LD(ALPH,ALFB,AH8,BH8,CH8,DH8,HCNAX,ANS) 
IF(SFAC(L).LE. 0.000)  GO  TO  62 
APLOG-LOG (SFAC( L) *C0N1*FAC*ANS*ZTT/TAU2 ( L ) ) 

62  FCORR-0.000 

IF(IFAN.Eq.2)  GO  TO  65 

DFAN-LOG( FCON*BETQ(L)/(ZFT* (1 . ODO-BETQ(L) ) ) )+UINF-0 . 5D0*BETQ( L) 
FCORR-DFAN/Z 
65  GRANO(L)-0.0 

IF(SFAC(L).LE. 0.000)  GO  TO  70 

GRANO( L ) -TAUFAC (L ) *SFAC ( L ) * ( APLOG-FCORR ) 

70  CONTINUE 

CALL  GRALD(OELTA, GRAND, LNAX,RES2) 

APL(H)-C0NZ«RES2 
80  CONTINUE 
CHICQ-0.0 
00  90  H-1,HHAX 

90  CHICQ-CHICQ+CP(H) 

CHIAQL-O.O 

00  91  M-1,HHAX 

91  CHIAQL-CHIAQL+APL(H) 

CHIAQL-CHIAQL/CHICQ 

CHIAQ-EXP(CHIAQL) 

COL-CHICQ/CHIAQ 

SB-LOG( COL ) +1 . 000-2 . ODO«EUL 
B -1 . 1500+ 1 . 12200*L0G ( COL ) 

110  BNEW-B-B«(B-LOG(B)-SB)/(B-1.0) 

ARG-ABS((BNEW-B)/B) 

IF(ARG-BEPS)130,120,120 
120  B-BNEW 
GO  TO  110 
130  B-BNEW 

SCALE-SQRT(CHICq«B) 

DHORM- ( SCALE*PATH ) **2 
RMAX-PATH*SCALE*THMAX 
DX-RNAX/DBLE(KNAX-1) 

OELTX-DX/3.0D0 
DO  240  K-1,KHAX 
X(K)-DX*REAL(K-1) 

TH-(X(K)/PATH) /SCALE 
IF(TH-THCUTl) 165, 165, 170 

165  CALL  BSP0L0(TH,THA1,AS1,BS1,CS1,DS1,JHA1,C0F1(K)) 

GO  TO  180 

170  CALL  BSPOLO(TH,THB1,AL1,BL1,CL1,OL1,JHB1,RES1) 

TERH-0.5D0*(TH*M) 

C0F1(K)-RES1/TERM 
180  IF(TH-THCUT2) 190, 190,200 

190  CALL  BSP0L0(TH,THA2,AS2,BS2,CS2,0S2,JHA2,C0F2(K)) 

GO  TO  210 

200  CALL  BSPOLO(TH,THB2,AL2,BL2,CL2,OL2,JMB2,RES2) 

TERH-0.500*(TH*M) 

COF2(K}-RES2/TERH 
210  ARG-TH*TH 

IF(ARG-86. 0)230, 230, 220 
220  COFO-0.000 
GO  TO  235 

230  C0F0-2.0D0*EXP(-ARG) 

235  DOR (K )- (COFO+COFl (K )/B+C0F2 (K ) /BB ) /DNORN 
0(K)-X(K)«D0R(K) 

240  DOR(K)-DOR(K)/(2.000«PI) 

CALL  Q<ALD(DELTX,D,KHAX,SUM) 

00  245  K-1,KMAX 
DOR(K)-OOR(K)/SUH 
245  0(K)-0(K)/SUN 

IF(N.NE.l)  GO  TO  257 
WRITE  (8,250) 

250  FORHAT('  TIN  - energjr  of  Incident  proton  bean,  HeV') 

WRITE  (8,251) 

251  FORHAT('  RGIN  - CSOA  range  at  energy  TIN,  g/cin2') 

WRITE  (8,252) 

252  FORHATC'  ZR  - depth.  In  units  of  RGIN') 

WRITE  (8,253) 

253  FORHAT('  RHAX  - largest  radial  distance') 

WRITE  (8,254) 

254  FORNAT('  SUM  - cumulative  Integral  of  Nollere  distribution  up  to 
IRKAX') 

WRITE  (8,1) 

WRITE  (8,255) 

255  F0RHAT(9X,' TIN', 8X, 'RGIN') 

WRITE  (8,256)  TIN, RGIN 

256  FORMAT(2F12.6) 

WRITE  (8,1) 


257  WRITE  (8,258) 

258  FORMAT(10X, 'ZR' ,8X, 'RMAX' ,9X, 'SUM' ) 

WRITE  (8,260)  ZR(N) ,RMAX,SUH 

260  FORMAT (3F12. 6) 

WRITE  (8,265) 

265  FORMAT( 'Radial  distances,  rho,  cm') 

WRITE  (8,270)  (X(K) ,K-1,KMAX) 

WRITE  (8,275) 

275  FORMAT ('Radial  distribution,  f(rho),  cm-2') 

WRITE  (8,270)  (DOR(K) ,K-1,KMAX) 

270  F0RMAT(1P8E12.5) 
print  *,n 
280  CONTINUE 
STOP 
END 

SUBROUTINE  SCOFD(X,F,NMAX,A,B,C,0) 

C 1 JAN  87.  DOUBLE  PRECISION 
C IF  S LIES  BETWEEN  X(M)  AND  X(M+1),  THEN 
C F(S)-( (D(M)*S+C(N) )*S+B(M) )*S+A(M) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  X(IOOO) ,F(1000) ,A(1000) ,B(1000) ,C(1000) ,0(1000) 
Nl-2 

M2-NNAX-1 

S-0.0 

DO  10  N-1,M2 
0(M)-X(M+1)-X(N) 

R-(F(M+1)-F(M))/D(M) 

C(N)-R-S 
10  S-R 
S-0.0 
R-0.0 
C(l}-0.0 
C(NNAX)-0.0 
DO  20  N-N1,M2 
C(H)-C(M)+R*C(M-1) 

B(M)-(X(M-1)-X(M+1))*2.0D0-R*S 

S-D(N) 

20  R-S/B(N) 

NR-N2 

DO  30  N-N1,M2 

C(MR)-(D(MR)*C(MR+1)-C(MR))/B(MR) 

30  MR-MR-l 
00  40  N-1,M2 
S-O(N) 

R-C(N+1)-C(N) 

0(N)-R/S 

C(N)-C(M)*3.0D0 

B(N)-(F(M+1)-F(M))/S-(C(M)+R)*S 
40  A(N)-F(N) 

RETURN 

END 

SUBROUTINE  BSPOLO(S,X,A,B,C,D,N,G) 

C 1 JAN  87.  DOUBLE  PRECISION 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  X(1000),A(1000),B(1000),C(1000), 0(1000) 

IF  (X(l).GT.X(N))  GO  TO  10 
IDIR-0 
MLB-0 
MUB-N 
GO  TO  20 
10  IDIR-1 
MLB-N 
MUB-O 

20  IF  (S.GE.X(NUB+IDIR))  GO  TO  60 
IF  (S.LE.X(NLB+1-IDIR))  GO  TO  70 
ML-MLB 
MU-MUB 
GO  TO  40 

30  IF  (lABS(MU-ML).LE.l)  GO  TO  80 
40  NAV-(NL+NU)/2 

IF  (S.LT.X(MAV))  GO  TO  50 
ML-HAV 
GO  TO  30 
50  MU-MAV 
GO  TO  30 

60  NU-NUB+2* IDIR-1 
GO  TO  90 

70  HU-NLB-2«IDIR+1 
GO  TO  90 
80  MU-MU+IOIR-l 
90  Q-S-X(NU) 

G-( (D(NU)*q+C(HU) )*q+B(MU) )*q+A(MU) 

RETURN 

END 

SUBROUTINE  GRAL0( DELTA, G,N, RESULT) 

C 1 JAN  87.  DOUBLE  PRECISION 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  G(2001) 

NLl-N-1 

NL2-N-2 

IF  (DBLE(N)-2.0D0*DBLE(N/2))  100,100,10 
C IF  N IS  ODD,  GO  TO  10  - IF  N IS  EVEN,  GO  TO  100 
10  IF(N-1)15, 15,20 
15  SIGMA-O.O 
GO  TO  70 

20  IF(N-3)30,30,40 


72 


30  SIGMA-G(1)+4.0D0*G(2)+G(3) 

GO  TO  70 
40  SUH4-0.0 

DO  50  K-2,NL1,2 
50  SUH4-SUH4>G(K) 

SUH2-0.0 
DO  60  K-3,NL2,2 
60  SUH2-SUH2>G(K) 

SIGHA-G(  1 )44. 0D0«SUH44^2 .0D0«SUM2-t^G(N) 

70  RESULT-DELTA*SIGHA 
RETURN 

100  IF(N-2)110,110,120 
110  SIGHA-1.500*(G(l)i-G(2)) 

GO  TO  70 

120  IF(N-4)130,130,140 

130  SIGHA-1.125DO*(G(1)+3.0DO*G(2)+3.0DO*G(3)+G(4)) 

GO  TO  70 

140  IF(N-6)150,150,160 

150  S IGMA-G ( 1 ) +3 . 875D0*G ( 2 ) +2 . 625DO*G ( 3 ) +2 . 625D0*G ( 4 ) + 
1 3.875D0*G(5)>G(6) 

GO  TO  70 

160  IF  (N-8)170,170,180 

170  S IGMA-G (1)^3. 875D0*G ( 2 ) +2 . 62500*G (3 ) *2 . 62500*G( 4) ^ 
1 3.875DO*G(5)+2.000*G(6)+4.000*G(7)+G(8) 

GO  TO  70 

180  SIG6-G( 1)^3. 875D0«G ( 2 ) +2 . 625D0«G( 3 ) *2 . 62500«G (4 ) ^ 

1 3.875D0*G(S)4^G(6) 

SUM4-0.0 

DO  190  K-7,NL1,2 
190  SUH4-SUH44G(K) 

SUH2-0.0 

DO  200  K-8,NL2,2 
200  SUH2-SUH24^G(K) 

S IGHA-S IG6-I-G  ( 6 ) >4 . 0D0*SUM4^2 . 000*SUH2  4^0  ( N ) 

GO  TO  70 
END 
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WRITE  (8,133) 

133  FORMAT('  TIM  - Energy  of  Incident  proton  beam,  MeV') 

WRITE  (8,134) 

134  FORMAT('  BR  - depth,  in  units  of  CSDA  range  at  energy  TIN') 
WRITE  (8,135) 

135  FORMAT('  RTOP  - maximum  radial  distance') 

WRITE  (8,1) 

136  WRITE  (8,138) 

138  FORMAT(9X, 'TIN' ,10X, 'BR' ,8X, 'RTOP' ) 

WRITE  (8,140)  TIN,BR,RTOP 
140  F0RMAT(4F12.6) 

DO  145  J-2,JHAX 
145  F(J)-F(J)/(RHO(J)*2.0*PI) 

F(1)-2.0*F(2)-F(3) 

WRITE  (8,148) 

148  F0RHAT(' Radial  distances,  rho,  cm') 

WRITE  (8,150)  (RHO(J),J-l,JMAX) 

WRITE  (8,149) 

149  FORHAT('Radiai  distribution,  f(rho),  cn-2') 

WRITE  (8,150)  (F(J),J-1,JHAX) 

150  F0RHAT(1P8E12.5) 

160  CONTINUE 

STOP 

END 


00  10  H-1,HHAX 
10  EBEGL(H)-LOG(EBEG(H)) 

OPEN  (UNIT-7, FILE-'STOPRANG. WAT') 

READ  (7,*)  KMAX 

READ  (7,«)  (ENG(K),K-1,KNAX) 

READ  (7,*)  (STP(K),K-1,KHAX) 

READ  (7,*)  (RANG(K),K-1,KKAX) 

CLOSE  (7) 

DO  20  K-1,KHAX 
EN6L(K)-LOG(ENG(K)) 

20  RAN6L(K}-L0G(RANG(K)4l. 00-09) 

CALL  SC0F(ENGL,RANGL,KHAX,A1,B1,C1,D1) 
TINL-LOG(TIN) 

CALL  BSP0L(TINL,ENGL,A1,B1,C1,D1,KHAX,RES) 
RGIN-EXP(RES) 

IF(TIN-EBEG(1)}30,30,40 
30  IF(TIN-EBE6(HHAX))40,60,60 
40  PRINT  SO 

50  F0RNAT(' Input  energy  is  out  of  range.') 

STOP 

60  DO  80  H-2,HH/U( 

IF(TINL-EBEGL(H) )80, 70,70 
70  Hl-M-1 
H2-N 
GO  TO  85 
80  CONTINUE 

85  FAC1-(TINL-EBEGL(H2))/(EBEGL(H1)-EBEGL(H2)) 
FAC2-1.0-FAC1 
DO  160  L-1,LHAX 
CAa  INOEX(LP(L},EXT) 

IWHJT-'RADF. '//EXT 
OPEN  (7, INPUT) 

READ  (7,5)  LINE 
READ  (7,*)  JHAX,TI1AX 
READ  (7,*)  BR 
DO  90  J-1,JHM( 

READ  (7,*)  IHD,T{J),(FIT(J,M),M-1,MHAX) 

90  CONTINUE 
CLOSE  (7) 

IF(L.6T.l)  GO  TO  95 
95  DELTA-T(JHAX)/REAL(JHAX-l)/3.0 
DO  110  H-1,HHAX 
DO  100  J-1,JNAX 
100  FT(J)-FIT(J,H) 

CALL  GRAL(DELTA,FT,JHAX,SUH(N)) 
no  CONTINUE 

DO  120  H-1,NHAX 
DO  120  J-1,JHAX 
120  FIT(J,N)-FIT(J,H)/SUN(H) 

DO  130  J-1,JHAX 

F(J)-{FAC1*FIT(J,M1)+FAC2*FIT(J,H2))/RGIN 

130  RHO(J)-RGIN*T(J) 

OELTA-OELTA*RGIN 

CAU  6RAL(DELTA,F,JHAX,SUNH) 

RTQP-THAX*RGIN 
IF(L.NE.l)  GO  TO  136 
WRITE  (8,131)  LMAX 

131  F0RMAT(I6,'  depths') 

WRITE  (8,132)  JHAX 

132  F0RHAT(I6,'  radial  distances') 

WRITE  (8,1) 


PROGRAM  PTRAD 

C 

C 20  Apr  93.  Written  by  Martin  J.  Berger,  NIST. 

C 

C Calculates  radial  dose  distribution  f(rho)  for  protons  in 

C water  at  depths  equal  to  the  following  fractions  of  the 

C CSDA  range:  0.985,  0.99,  0.995,  1.0,  1.005,  1.01, 

C 0.015,  1.020,  1.025,  for  beam  energies  between 

C 50  and  250  HeV.  Results  are  obtained  by  interpolation 

C in  a database  generated  with  the  Monte  Carlo  program 

C PTRAfOD. 

C 

C The  following  input  files  must  be  available: 

C STOPRANG.WAT:  Stopping-power  and  range  table  for  water. 

C RAOF.032,..,RAOF.04O:  database  of  processed  Monte 

C Carlo  results  from  PTRAN3D. 

C 

DIMENSION  EBEG(7),EBEGL(7),T(401),RH0(401),FIT(401,7), 

1 F(401),ENG(133),ENGL(133),STP(133),RAN6(133),RANGL(133), 

2 A1(133),B1(133),C1(133),D1(133),FT(401),SUM(7),LP(9) 

CHARACTER  EXT*3,INPUT*16,OUTPUT*16 

DATA  LP/32,33,34,35,36,37,38,39,40/,MMAX/7/,LHAX/9/ 

DATA  EBEG/250. 0,200. 0,160. 0,130. 0,100. 0,70. 0,50.0/ 

1 F0RMAT(1H  ) 

5 FORHAT(A) 

PI-4.0*ATAN(1.0) 

PRINT  *,'  Enter  energy  (MeV):  ' 

READ  *,  TIN 

PRINT  *,'  Suggested  name:  PTRAD. nnn' 

PRINT  Enter  name  of  output  file:  ' 

READ  5,  OUTPUT 
OPEN  (8,0UTPUT) 

WRITE  (8,6)  OUTPUT 

6 FORMAT ('Program  PTRAD,  output  file  ',A) 

WRITE  (8,1 
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PROGRAH  FCIR 


C 

C 15  May  93.  Written  by  Martin  J.  Berger,  NIST. 

C 

C Calculates  reduction  factors  and  absorbed-dose  values 

C for  a circular  field,  as  a function  of  the  radial 

C distance  from  the  center  of  the  field. 

C 

C Calls  subroutines  CIRCLE,  SCOFO,BSPOLD  and  6RAL0. 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  R(401),F(401),FL(401),ACOF(401),BCOF(401), 

1 CCOF (401 ) , DCOF (401 ) , GRAN0( 5121 ) , S ( 401 ) , REDUC (401 ) , DRAO(401 ) , 

2 ZR(120),OODX(120),REO(120),OAX(120), 

3 ZD(100),DEPOT(100) ,AP(100) ,BP(100) ,CP(100) ,0P(100) 

CHARACTER*16  INPUTl, INPUT2, INPUT3, OUTPUT, INFIL 
CHARACTER  LINE«80 

1 F0RMAT(1H  ) 

DATA  NNX/161/,EREDUCL/1.0D-04/ 

5 FORMAT(A) 

P 1-4. 0D0*AT AN (1.000) 

PRINT  Enter  radius  of  circular  field  (cm):  ' 

READ  *,  RADIUS 

PRINT  Options  for  specifying  radial  distances' 

PRINT  1)  Get  set  of  values  from  file' 

PRINT  2)  Get  set  of  values  from  keyboard' 

PRINT  *,'  3)  Specify  In  terns  of  maximum  distance' 

PRINT  and  number  of  distances:  ' 

PRINT  Enter  choice:  ' 

READ  INS 

GO  TO  (6,7,8),  INS 

6 PRINT  *,'  Enter  name  of  file  with  distance  values:  ' 

READ  5,  INFIL 

OPEN  (7, INFIL) 

READ  (7,*)  LMAX 

READ  (7,*)  (S(L),L-1,LMAX) 

CLOSE  (7) 

GO  TO  9 

7 PRINT  ♦,'  Enter  number  of  radial  distances  (no  greater  than  401): 

1' 

READ  *,  LMAX 

PRINT  Enter  distances  (cm):  ' 

READ  *,  (S(L),L-1,LMAX) 

GO  TO  9 

8 PRINT  Enter  maximum  radial  distance  (cm):  ' 

READ  *,  SHAX 

LHAX-1 

IF(SNAX.LE.O.ODO)  GO  TO  9 

PRINT  Enter  number  of  radial  distances  (no  greater  than  401): 
1' 

READ  «,  LMAX 

9 PRINT  Enter  Input  file  1 (from  MORAD):  ' 

READ  5,  INPUTl 

PRINT  *,'  Options  for  depths:' 

PRINT  *,'  1}  No  depths  greater  than  981  of  CSDA  range' 

PRINT  2)  Some  depths  greater  than  98t  of  CSDA  range' 

PRINT  Enter  choice:  ' 

READ  *,  INDEP 
IF(INDEP.EQ,1)  GO  TO  10 

PRINT  *,'  Enter  Input  file  2 (from  PTRAD):  ' 

READ  5,  INPUT2 

10  PRINT  *,'  Enter  Input  file  3 (from  PTPOL):  ' 

READ  5,  INPUTS 

PRINT  *,'  Enter  name  of  output  file:  ' 

READ  5,  OUTPUT 
IF(INS.LT.3)  GO  TO  20 
S(l)-0.0 

IF(LMAX.LE.l)  GO  TO  20 
DS-SNAX/0BLE(LMAX-1) 

DO  IS  L-2,LHAX 
IS  S(L)-DS*0BLE(L-1) 

20  OPEN  (8,0UTPUT) 

WRITE  (8,25)  OUTPUT 

25  FORMAT ('Program  FCIR,  output  file  ',A) 

WRITE  (8,26)  INPUTl 
IF(INDEP.EQ.2)  WRITE  (8,26)  INPUT2 
WRITE  (8,26)  INPUTS 

26  FORMAT('Input  file  ',A) 

IF( INDEP. EQ.l)  WRITE  (8,1) 

WRITE  (8,1) 

WRITE  (8,27) 

27  F0RMAT('  TIN  • energy  of  Incident  proton  beam,  MeV') 

WRITE  (8,28) 

28  F0RKAT('  RGIN  - CSDA  range  at  energy  TIN,  g/cm2') 

WRITE  (8,30) 

30  F0RMAT('  RADIUS  • radius  of  circular  field,  cm') 

WRITE  (8,31) 

31  F0RMAT('  LMAX  • number  of  radial  distances  from  center  of  field' 

1) 

WRITE  (8,32) 

32  F0RNAT('’  NCASE  - number  of  depths') 

WRITE  (8,1) 

OPEN  (6, INPUTl) 

IF(INDEP.EQ.2)  OPEN  (7,INPUT2) 

READ  (6,5)  LINE 
READ  (6,5)  LINE 


READ  (6,*)  NCASEl 
READ  (6,*)  JMAXl 
DO  33  LN-1,8 
READ  (6,5)  LINE 

33  CONTINUE 

READ  (6,*)  TIN, RGIN 
READ  (6,5)  LINE 
NCASE2-0 
JHAX2-0 

IF( INDEP. EQ.l)  GO  TO  35 
READ  (7,5)  LINE 
READ  (7,5)  LINE 
READ  (7,*)  NCASE2 
READ  (7,*)  JMAX2 
DO  34  LN-1,5 
READ  (7,5)  LINE 

34  CONTINUE 

35  NCASE-NCASEl4^NCASE2 
WRITE  (8,36) 

36  F0RMAT(9X,' TIN', 8X, 'RGIN', 6X, 'RADIUS','  LMAX  NCASE') 
WRITE  (8,37)  TIN, RGIN, RADIUS, LMAX, NCASE 

37  F0RMAT(3F12. 5,216) 

WRITE  (8,1) 

IF(LMAX.NE.l)  WRITE  (8,38) 

38  F0RMAT( 'Radial  distances  from  center  of  field  (cm)') 
IF(LMAX.NE.l)  WRITE  (8,39)  (S(L),L-1,LMAX) 

39  F0RHAT(8F12.6) 

0PEN(11, INPUTS) 

DO  40  LN-1,9 
READ  (11,5)  LINE 

40  CONTINUE 

READ  (11,*)  G1,G2,G3,G4,KMAX 
READ  (11,5)  LINE 
READ  (11,5)  LINE 
READ  (11,*)  (ZD(K),K-1,KHAX) 

READ  (11,*)  LINE 

READ  (11,*)  (DEPOT(K),K-l,KMAX) 

CLOSE  (11) 

CALL  SCOFD(ZD,DEPOT,KHAX,AP,BP,CP,DP) 

DO  130  NC-1, NCASE 
IF(NC-NCASE1)42,42,43 

42  JHAX-JMAXl 
READ  (6,5)  LINE 

READ  (6,*)  ZR(NC),RMAX,SUH 

READ  (6,5)  LINE 

READ  (6,*)  (R(J),J-1,JMAX) 

READ  (6,5)  LINE 

READ  (6,*)  (F(J),J-1,JMAX) 

GO  TO  46 

43  JHAX-JMAX2 
READ  (7,5)  LINE 

READ  (7,*)  TIN,ZR(NC),RMAX 

READ  (7,5)  LINE 

READ  (7,*)  (R(J),J-1,JMAX) 

READ  (7,5)  LINE 

READ  (7,*)  (F(J),J-1,JMAX) 

JTAL-0 

DO  44  J-1,JMAX 
IF(F(J))45,45,44 

44  JTAL-JTAL+1 

45  JHAX-JTAL 

46  DO  47  J-1,JMAX 

47  FL(J)-LOG(2.0D0*PI*F(J)) 

CALL  SCOFD( R , FL , JMAX , ACOF , B COF , CCOF , DCOF ) 

DO  90  L-1,LHAX 

NMAX-NMX 

OLOVAL-2.000 

48  IF(S(L).GE. RADIUS)  GO  TO  70 
RHHAX-RAOIUS-S(L) 

I F ( RMAX-RHHAX ) 49 , 49 , 50 

49  REDUC(L)-1.0 
GO  TO  90 

50  DRH-RKHAX/DBLE(NHAX-1) 

DO  51  N-1,NHAX 
RHO-ORH*DBLE(N-l) 

CALL  BSPOLO ( RHO , R , ACOF , BCOF , CCOF , DCOF , JMAX, RES ) 

51  6RANO(N)-RHO*EXP(RES) 

DELTA-ORH/3.0D0 

CALL  6RALO(DELTA,GRAND,NHAX, PARTI) 

IF(S(L))52,52,S3 

52  REDUC(L) -PARTI 
GO  TO  90 

53  RHNIN-RAOIUS-S(L) 

RHHAX-RA0IUS4-S(L) 

RHNAX-N IN ( RHMAX , RHAX ) 

DRH-(RHHAX-RHHIN)/DBLE(NMAX-1 ) 

DO  55  N-1,NHAX 
RHO-RmiN4-DRH*OBLE(N-l ) 

CALL  BSPOLO(RHO,R, ACOF, BCOF, CCOF,DCOF,JMAX,RES) 
FINT-RHO*EXP(RES) 

CALL  CIRCLE(S(L), RADIUS, RHO, FRAC) 

55  GRANO(N)-FINT*FRAC 
DELTA-DRH/3.0D0 

CALL  6RALD(0ELTA,GRAND,NMAX,PART2) 

REDUC(L)-PART1+PART2 

GO  TO  90 
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70  RHHIN-S(L)-RAOIUS 
RHHAX-S(L}-f  RADIUS 
RHHAX -H IN ( RHHAX , RHAX ) 

DRH-(RHHAX-RHHIN)/DBLE(NHAX-1 ) 

DO  80  N-1,NHAX 
RHO-RHH  IN4^DRH*DBLE  ( N- 1 ) 

CALL  BSPOLD ( RHO , R , ACOF , BCOF , CCOF , DCOF , JMAX , RES ) 
FINT-RHO*EXP(RES) 

CALL  CIRCLE(S(L), RADIUS, RHO, FRAC) 

80  GRAND(N}-FINT*FRAC 
0ELTA-DRH/3.0D0 

CALL  GRALO(DELTA, GRAND, NHAX,REDUC(L)) 
DIFF-ABS(OLDVAL-REDUC(L} ) 

IF(DIFF>EREDUCL)90, 81,81 

81  OLDVAL-REDUC(L) 

IF(NHAX-2SS1)82,82,90 

82  NNAX-2*NHM(-1 
GO  TO  48 

90  CONTINUE 

CALL  BSPOLO(ZR(NC},ZD,AP,BP,CP,DP,KHAX,DODX(NC)) 

DO  95  L-1,LHAX 
95  DRAD(L)-DDDX(NC)«REDUC(L) 

IF(LMAX.NE.l)  WRITE  (8,100)  ZR(NC) 

100  FORHAT('Depth  as  fraction  of  RGIN  - ',F5.3) 

IF(LMAX.Eq.l)  REO(NC)-REDUC(l) 

IF(LNAX.EQ.l)  0AX(NC)-DRAD(1) 

IF(LHAX.NE.l)  WRITE  (8,105) 

105  FORMAT ('Reduction  factors') 

IF(LMAX.NE.l)  WRITE  (8,110)  (REDUC(L) ,L-1,LMAX) 

110  FORMAT(8F12.6) 

IF(LMAX.NE.l)  WRITE  (8,115) 

115  FORMAT ( 'Absorbed  dose,  MeV/g,  from  a beam  of  1 proton/cm2') 
IF(LMAX.NE.l)  WRITE  (8,110)  (DRAD(L),L-1,LHAX) 
print  120,  ncase,nc 

120  format('  ncase  - ',13,'  nc  - ',13) 

130  CONTINUE 

IF(LNAX.NE.l)  GO  TO  170 
WRITE  (8,140)  S(l) 

140  FORMAT ('Radial  distance  (cm)  from  center  of  field  - ',F5.3) 
WRITE  (8,1) 

WRITE  (8,150) 

150  FORMAT ('Depths,  In  units  of  RGIN') 

WRITE  (8,160)  (ZR(NC),NC-1, NCASE) 

160  FORNAT(8F12.6) 

WRITE  (8,105) 

WRITE  (8,160)  (RED(NC),NC-1, NCASE) 

WRITE  (8,115) 

WRITE  (8,160)  (DAX(NC),NC-1, NCASE) 

170  STOP 
END 

SUBROUTINE  CIRCLE (S, RADIUS, RHO, FRAC) 

C 4 Feb  93.  Calculates  arclength  fraction  In  circle. 

C 

C Circle  has  center  at  origin  and  has  radius  R. 

C Reference  point  Is  at  x>-0  ,ir-0. 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

PI-4.0D0«ATAN(1.000) 

IF(S-RAOIUS)20, 10,20 
10  ARG-RHO/(2.000*RAOIUS) 

GO  TO  30 

20  ARG-( S**2+RHO**2-RAOIUS**2 ) / ( 2 . ODO*S*RHO) 

30  IF(ARG.LT.-l.ODO)  ARG— l.ODO 
IF(ARG.GT.l.ODO)  ARG-l.ODO 
FRAC-ACOS(ARG)/PI 
RETURN 
END 
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CALF-COS(ALF) 

SALF-SIN(ALF) 

DO  205  L-1,LHAX 
NHAX-NHX 
OLDVAL-2.000 
X-CALF*DS*DBLE(L-1) 

Y-SALF*DS*DBLE(L-1) 

100  IF(X.LE.A.AND.Y.LE.B)  GO  TO  110 
IF(X.LE.A.AND.Y.GT.B)  GO  TO  140 
IF(X.GT.A.AND.Y.LE.B)  GO  TO  160 
IF{X.GT.A.AND.Y.GT.B)  GO  TO  180 
110  RHHAX-HIN(A-X,B-Y) 

IF(RMAX-RHMAX)115,115,120 
115  REDUC(L) -1.000 
GO  TO  200 

120  ORH-RHHAX/DBLE(NHAX-l) 

DO  125  N-1,NHAX 
RHO-ORH*OBLE(N-l) 

CALL  BSPOLD ( RHO , R , ACOF , BCOF , CCOF , DCOF , JHAX , RES } 
125  GRAND(N)-RHO*EXP(RES) 

DELTA-DRH/3.000 

CALL  GRAL0( DELTA, GRAND, NHAX, PARTI) 
RHHIN-HIN(A-X,B-Y) 

RHMAX-SQRT((X+A)**2+(Y+B)**2) 

RHMAX-HIN(RHHAX,RHAX) 

0RH-(RHHAX-RHHIN)/DBLE(NHAX-1) 

DO  130  N-1,NHAX 
RHO-RHHIN«ORH«DBLE( N- 1 ) 

CALL  BSPOLD ( RHO, R , ACOF , BCOF , CCOF , DCOF , JHAX , RES ) 
FINT-RHO*EXP(RES) 

CALL  CORNER(RHO,A^X,B+Y,FRAC1) 

CALL  CORNER ( RHO, A4X,B-Y,FRAC2) 

CALL  C0RNER(RH0,A-X,BW,FRAC3) 

CALL  C0RNER(RH0,A-X,B-Y,FRAC4) 

1 30  GRAND ( N ) -F INT* ( FRACl + FRAC2+ FRAC3+ FRAC4 ) 
DELTA-DRH/3.0D0 

CALL  GRALD( DELTA, GRAND, NHAX, PART2) 
REDUC(L)-PARTl4PART2 
GO  TO  200 
140  RHHIN-Y-B 

IF(RHAX-RHHIN)142,142,145 
142  REDUC(L)-0.0D0 
GO  TO  200 

145  RHMAX-SQRT((X+A)**2+(Y+B)**2) 

RHHAX-H IN ( RHNAX , RHAX ) 
DRH-(RHHAX-RHHIN)/0BLE(NHAX-1) 

DO  150  N-1,NHAX 
RH0-RHNIN+0RH«DBLE(N-1) 

CALL  BSPOLD ( RHO , R , ACOF , BCOF , CCOF , DCOF , JHAX , RES ) 
FINT-RHO*EXP(RES) 

CALL  CORNER(RHO,A+X,Y+B,FRACl) 

CALL  C0RNER(RH0,A-X,Y^B.FRAC2) 

CALL  C0RNER(RH0,A+X,Y-B,FRAC3} 

CALL  C0RNER(RH0,A-X,Y-B,FRAC4) 

150  GRAH0(N)-FINT*(FRAC1+FRAC2-FRAC3-FRAC4) 
0ELTA-DRH/3.0D0 

CALL  GRALD( DELTA, GRAND, NHAX,REDUC(L)) 

GO  TO  200 
160  RHHIN-X-A 

I F ( RHAX -RHH IN ) 1 62 , 1 62 , 1 65 
162  REDUC(L)-0.0D0 
GO  TO  200 

165  RHHAX-SQRT((X+A)**2+(Y+B)**2) 
RHHAX-HIN(RHHAX,RHAX) 
DRH-(RHNAX-RHHIN)/DBLE(NMAX-1 ) 

DO  170  N-1,NNAX 
RHO-RHN IN^DRH*OBLE ( N- 1 ) 

CALL  BSPOLD(RHO,R, ACOF, BCOF, CCOF, DCOF, JHAX,RES) 
FINT-RHO*EXP(RES) 

CALL  CORNER(RHO,X^A,B+Y,FRACl) 

CALL  CORNER(RHO,X+A,B-Y,FRAC2) 

CALL  C0RNER(RH0,X-A,B*Y,FRAC3) 

CALL  C0RNER(RH0,X-A,B-Y,FRAC4) 

170  GRAN0(H)-FIHT*(FRAC1+FRAC2-FRAC3-FRAC4) 
DELTA-DRH/3.000 

CALL  6RALD(0ELTA, GRAND, NHAX, REOUC(L)} 

GO  TO  200 

180  RHHIM-SQRT((X-A)**2+{Y-B)**2) 

IF(RHAX-RHHIN)182, 182,185 
182  REDUC(L) -0.000 
GO  TO  200 

185  RHHAX-SQRT((X+A)**2+(Y+B)**2) 

RIWAX-H  IN  ( RmAX , RHAX ) 

DRH-<RHHAX-RHHIN) /DBLE(NHAX-1 } 

DO  190  N-1,NHAX 
RH0-RHHIN+0RH*0BLE(N-1) 

CALL  BSPOLO(RHO,R, ACOF, BCOF, CCOF, DCOF, JHAX, RES) 
FINT-RHO«EXP(RES) 

CALL  CORNER ( RHO, X4^A,Y4B,FRAC1) 

CALL  C0RNER(RH0,X-A,Y4^B,FRAC2) 

CALL  C0RNER(RH0,X>A,Y-B,FRAC3) 

CALL  C0RNER(RH0,X-A,Y-B,FRAC4) 

190  GRAN0(N)-FINT*(FRAC1-FRAC2-FRAC34FRAC4) 
DELTA-0RH/3.0D0 

CALL  GRAL0( DELTA, GRAND, NHAX, REOUC(L)) 

200  OIFF-ABS(OLOVAL-REOUC(L)) 


IF(DIFF-EREDUCL)205,201,201 

201  OLDVAL-REOUC(L) 

I F ( NHAX-2561 ) 202 , 202 , 205 

202  NHAX-2«NHAX-1 
GO  TO  100 

205  CONTINUE 

CALL  BSPOLO(ZR(NC),ZD,AP,BP,CP,OP,KHAX,DOOX(NC)) 

DO  206  L-1,LHAX 

206  DRAO(L)-DDOX(NC)«REDUC(L) 

IF(LHAX.NE.l)  WRITE  (8,210)  ZR{NC) 

210  FORHAT('Oepth  as  fraction  of  RGIN  - ',F5.3) 

IF(LHAX.EQ.l)  RED(NC)-REDUC(1) 

IF(LHAX.EQ.l)  0AX(NC)-DRA0(1) 

IF(LHAX.NE.l)  WRITE  (8,220) 

220  FORHAT( 'Reduction  factors') 

IF(LHAX.NE.l)  WRITE  (8,230)  (REDUC(L) ,L-1,LHAX) 

230  F0RHAT(8F12.6) 

IF(LHAX.NE.l)  WRITE  (8,240) 

240  FORHAT( 'Absorbed  dose,  HeV/g,  from  a beam  of  1 proton/cm2') 
IF(LHAX.NE.l)  WRITE  (8,230)  (DRAO(L),L-l,LHAX) 
print  245,  ncase,nc 

245  format('  ncase  - ',13,'  nc  - ',13) 

250  CONTINUE 

IF(LHAX.NE.l)  GO  TO  290 
WRITE  (8,260)  S(l) 

260  FORHAT( 'Radial  distance  (cm)  from  center  of  field  - ',F5,3) 

WRITE  (8,1) 

WRITE  (8,270) 

270  FORHAT( 'Depths,  In  units  of  RGIN') 

WRITE  (8,280)  (ZR(NC) ,NC-1, NCASE) 

280  F0RHAT(8F12.6) 

WRITE  (8,220) 

WRITE  (8,280)  (RED(NC) ,NC.l, NCASE) 

WRITE  (8,240) 

WRITE  (8,280)  (DAX(NC),NC-1, NCASE) 

290  STOP 
END 

SUBROUTINE  CORNER(R,SA,SB,FRAC) 

C 27  Har  88.  Computes  rel .arc length  (fracton  of  2 pi)  for  point  In 
C corner  of  rectangle  with  sides  SA  and  SB. 

IHPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

TPI-8.000*ATAN( 1.000) 

IF(SA)60,60,5 

5 IF(SB)60,60,6 

6 A-HAX(SA,SB) 

B-HIN(SA,SB) 

IF(R-B)10,10,20 

10  FRAC-0.25DO 
RETURN 

20  IF(R-A)30,30,40 
30  FRAC-0.2500-ACOS(B/R)/TPI 
RETURN 

40  IF(R*R-A*A-8*B)50,50,60 
50  FRAC-0 . 25D0- ( ACOS ( B/R ) + ACOS ( A/R ) ) /TP I 
RETURN 

60  FRAC-0. 000 
RETURN 
END 
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PROGRAM  FREC 
C 

C 15  May  93.  Written  by  Martin  J.  Berger,  NIST. 

C 

C Calculates  reduction  factors  and  absorbed-dose  values 

C for  a rectangular  field. 

c 

C Center  of  coordinate  system  Is  at  x-0,  y-0. 

C Rectangle  extends  from  -A  to  A In  x,  -B  to  B In  y. 

c 

C Results  are  calculated  as  functions  of  the  distance 

C from  the  center  of  the  field,  along  a line  that  makes 

C an  angle  alpha  (between  0 and  90  degrees)  with  respect 

C to  the  x-axIs. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  REDUC(401),S(401),R(401),F(401),FL(401), 

1 ACOF ( 401 ) , BCOF ( 401 ) , CCOF ( 401 ) , DCOF ( 401 ) , GRAN0( 5121 ) , 

2 ZO(100),OEPOT(100),AP(100),BP(100},CP(100),DP(100),DRAD(401), 

3 ZR(401),0ODX(401),RED(4Ol),0AX(401) 

CHARACTER*16  INPUT1,INPUT2, INPUTS, OUTPUT, INFIL 
CHARACTER  LINE«80 

DATA  NMX/161/,EREDUCL/1.0D-04/ 

1 FORMAT (IH  ) 

5 FORMAT(A) 

PI-4.0DO*ATAN(1.0DO) 

RAOFAC-PI/180.0D0 

PRINT  Rectangular  field  extends  from  -A  to  A In  x,  -B  to  B In 
ly' 

PRINT  Enter  values  of  A and  6 (cm):  ' 

READ  «,  A,B 

PRINT  *,  ' Line  goes  through  a corner  of  the  rectangle  (1-yes,  2-n 
10):  ' 

READ  *,  JVEC 
IF(JVEC.NE.2)  GO  TO  10 

PRINT  Enter  angle  alpha  with  respect  to  x-axIs  (degrees):  ' 
READ  ♦,  ALPHA 
ALF-ALPHA*RAOFAC 
GO  TO  11 

10  ALF-ATAN(B/A) 

ALPHA-ALF/RAOFAC 

11  PRINT  Options  for  specifying  radial  distances:' 

PRINT  1)  Get  set  of  values  from  file' 

PRINT  *,'  2)  Get  set  of  values  from  keyboard' 

PRINT  3)  Specify  set  of  values  In  terms  of  maximum  distance' 

PRINT  *,'  and  number  of  distances,  LHAX:  ' 

PRINT  *,'  Enter  choice:  ' 

READ  *,  INS 

GO  TO  (12,13,14),  INS 

12  PRINT  *,'  Enter  name  of  file  with  distance  values:  ' 

READ  5,  INFIL 

OPEN  (7, INFIL) 

READ  (7,«)  LMAX 

READ  (7,*)  (S(L),L-1,LMAX) 

CLOSE  (7) 

GO  TO  15 

13  PRINT  *,'  Enter  ntnber  of  radial  distances  (no  greater  than  401): 
1' 

READ  *,  LHAX 

PRINT  *,'  Enter  distances  (cm):  ' 

READ  «,  (S(L),L-1,LHAX) 

GO  TO  15 

14  PRINT  Enter  maximum  radial  distance  (cm):  ' 

READ  «,  SMAX 

LHAX-1 

IF(SNAX.LE.O.ODO)  GO  TO  15 

PRINT  *,'  Enter  number  of  radial  distances  (no  greater  than  401): 
1' 

READ  *,  LHAX 

15  PRINT  *,'  Enter  name  of  input  file  1 (from  MORAD):  ' 

READ  5,  INPUTl 

PRINT  Options  for  depths:' 

PRINT  *,'  1)  No  depths  greater  than  98X  of  CSOA  range' 

PRINT  2)  Some  depths  greater  than  98S  of  CSDA  range' 

PRINT  *,'  Enter  choice:  ' 

READ  *,  INDEP 
IF(INDEP.EQ.l)  GO  TO  16 

PRINT  Enter  name  of  Input  file  2 (from  PTRAO):  ' 

READ  5,  INPUT2 

16  PRINT  *,'  Enter  name  of  Input  file  3 (from  PTPOL):  ' 

READ  5,  INPUT3 

PRINT  Enter  name  of  output  file:  ' 

READ  5,  OUTPUT 
IF(INS.LT.3)  GO  TO  18 
S(l)-0.0 

IF(LMAX.LE.l)  GO  TO  18 
DS-SNAX/DBLE( LHAX-1) 

DO  17  L-2,LHAX 

17  S(L)-0S*0BLE(L-1) 

18  OPEN  (8,0UTPUT) 

WRITE  (8,19)  OUTPUT 

19  F0RHAT( 'Program  FREC,  output  file  ',A) 

WRITE  (8, 20) INPUTl 
IF(INDEP.Eq.2)  WRITE  (8,20)  INPUT2 
WRITE  (8,20)  INPUT3 

20  F0RMAT( ' Input  file  - ',A) 


IF( INDEP. EQ.l)  WRITE  (8,1) 

WRITE  (8,1) 

WRITE  (8,25) 

25  F0RHAT('  TIN  - energy  of  Incident  proton  beam,  MeV') 

WRITE  (8,26) 

26  F0RMAT('  RGIN  - CSDA  range  at  energy  TIN,  g/cm2') 

WRITE  (8,27) 

27  FORMAT('  2A  - side  of  rectangular  field,  cm') 

WRITE  (8,28) 

28  FORMAT('  2B  - side  of  rectanglular  field,  cm') 

WRITE  (8,29) 

29  FORHAT('  ALPHA  • angle  (deg)  defined  below') 

WRITE  (8,30) 

30  F0RHAT('  LHAX  - number  of  radial  distances  from  center  of  field' 

1) 

WRITE  (8,31) 

31  F0RMAT('  NCASE  - number  of  depths') 

WRITE  (8,1) 

WRITE  (8,32) 

32  FORMAT ('Rectangular  field  extends  from  -A  to  A In  x,  and  from  -B  t 
lo  B In  y. ' ) 

WRITE  (8,33) 

33  F0RHAT('L1ne  along  which  reduction  factor  Is  calculated  starts  at 
Ithe  origin' ) 

WRITE  (8,34) 

34  FORHAT('and  makes  angle  ALPHA  (deg)  with  respect  to  x-ax1s.') 

WRITE  (8,1) 

OPEN  (6, INPUTl) 

IF(INDEP.EQ.2)  OPEN  (7,INPUT2) 

READ  (6,5)  LINE 
READ  (6,5)  LINE 
READ  (6,«)  NCASEl 
READ  (6,*)  JHAXl 
DO  35  LN-1,8 
READ  (6,5)  LINE 

35  CONTINUE 

READ  (6,*)  TIN, RGIN 
READ  (6,5)  LINE 
NCASE2-0 
JMAX2-0 

IF( INDEP. EQ.l)  GO  TO  37 
READ  (7,5)  LINE 
READ  (7,5)  LINE 
READ  (7,*)  NCASE2 
READ  (7,*)  JMAX2 
DO  36  LN-1,5 
READ  (7,5)  LINE 

36  CONTINUE 

37  NCASE-NCASE1+NCASE2 
WRITE  (8,38) 

38  F0RMAT(9X,'TIN',8X,'RGIN',11X,'A',11X,'B',7X,'ALPHA','  LMAX  NCASE 
1') 

WRITE  (8,39)  TIN,RGIN,A,B,ALPHA,LHAX,NCASE 

39  F0RMAT(5F12. 5,216) 

IF(LNAX.NE.l)  WRITE  (8,1) 

IF(LMAX.NE.l)  WRITE  (8,40) 

40  F0RMAT( 'Radial  distances  from  center  of  field  (cm)') 

IF(LMAX.NE.l)  WRITE  (8,41)  (S(L),L-1,LHAX) 

41  FORMAT(6F12.5) 

0PEN(11,INPUT3) 

DO  42  LN-1,9 
READ  (11,5)  LINE 

42  CONTINUE 

READ  (11,*)  61,62,G3,G4,KHAX 
READ  (11,5)  LINE 
READ  (11,5)  LINE 
READ  (11,*)  (2D(K),K-1,KMAX) 

READ  (11,*)  LINE 

READ  (11,*)  (DEPOT(K),K-l,KMAX) 

CLOSE  (11) 

CALL  SCOFD(ZD,DEPOT,KHAX,AP,BP,CP,DP) 

DO  250  NC-1, NCASE 
IF(NC-NCASE1)43,43,44 

43  JMAX-JHAXl 
READ  (6,5)  LINE 

READ  (6,*)  ZR(NC),RHAX,SUH 

READ  (6,5)  LINE 

READ  (6,*)  (R(J),J-1,JMAX) 

READ  (6,5)  LINE 

READ  (6,*)  (F(J),J-1,JHAX) 

GO  TO  47 

44  JMAX-JMAX2 
READ  (7,5)  LINE 

READ  (7,*)  TIN,ZR(NC),RMAX 

READ  (7,5)  LINE 

READ  (7,*)  (R(J),J-1,JMAX) 

READ  (7,5)  LINE 

READ  (7,*)  (F(J),J-1,JMAX) 

JTAL-0 

DO  45  J-1,JHAX 
IF(F(J))46,46,45 

45  JTAL-JTAL^l 

46  JHAX-JTAL 

47  DO  48  J-1,JHAX 

48  FL(J)-L0G(2.0D0*PI*F(J)) 

CALL  SC0FD( R, FL , JMAX, ACOF , BCOF , CCOF , DCOF ) 
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PROGRAM  FRECl 


15  Hay  93.  Written  by  Martin  J.  Berger,  NIST. 

Calculates  reduction  factors  and  absorbed-dose  values 
for  a rectangular  field. 

Center  of  coordinate  system  Is  at  x«0,  y»0. 

Rectangle  extends  from  -A  to  A In  x,  -B  to  B In  y. 

Results  are  calculated  a specified  set  of  (x,y)  values. 
Without  loss  of  generality,  It  Is  assumed  that  x and  y 
are  non-negative. 


IMPLICIT  DOUBLE  PRECISION  (A-H,0-2) 

DIMENSION  REDUC(401),XP(401),YP(401),R(401),F(401),FL(401), 

1 ACOF (401 ) , BCOF (401 ) , CCOF (401 ) , DCOF (401 ) , GRAND( 5121 ) , 

2 ZD(100),DEPOT(100),AP(100),BP(100),CP(100),DP(100),DRAD(401), 

3 ZR(401),DDDX(401) 

CHARACTER*16  INPUTl , INPUT2 , INPUT3 , OUTPUT , INFI L 

CHARACTER  LINE*80 

DATA  NHX/161/,EREDUCL/1.0D-04/ 

1 FORMAT (IH  ) 

5 FORHAT{A) 

PI-4.ODO«ATAN(1.0DO) 

PRINT  Rectangular  field  extends  from  -A  to  A In  x,  -B  to  B In 
ly' 

PRINT  Enter  values  of  A and  B (cm):  ' 

READ  *,  A,B 

PRINT  Options  for  specifying  (x,y)-values: ' 

PRINT  1)  Enter  values  from  keyboard' 

PRINT  2)  Enter  values  from  file:' 

PRINT  *,'  Enter  choice:  ' 

READ  «,  INP 
IF(INP.EQ.2)  GO  TO  7 

PRINT  *,'  Enter  number  of  (x,y)  values:  ' 

READ  *,  LMAX 
DO  6 L-1,LHAX 

PRINT  *,'  Enter  x and  y (cm,  non-negative):  ' 

READ  *,  XP(L),YP(L) 

6 CONTINUE 
GO  TO  9 


7 PRINT  *,'  Enter  name  of  input  file:  ' 

READ  5,  INFIL 

OPEN  (7, INFIL) 

READ  (7,*)  LMAX 

DO  8 L-1,LHAX 

READ  (7,*)  XP(L),YP(L) 

8 CONTINUE 
CLOSE (7) 

9 PRINT  *,'  Enter  name  of  Input  file  1 (from  MORAD) : ' 

READ  5,  INPUTl 

PRINT  *,'  Options  for  depths:' 

PRINT  *,'  1)  No  depths  greater  than  98*  of  CSDA  range' 

PRINT  *,'  2)  Some  depths  greater  than  98*  of  CSDA  range' 

PRINT  *,'  Enter  choice:  ' 

READ  *,  INDEP 
IF(INDEP.EQ.l)  GO  TO  16 

PRINT  *,'  Enter  name  of  Input  file  2 (from  PTRAD) : ' 

READ  5,  INPUT2 

16  PRINT  *,'  Enter  name  of  Input  file  3 (from  PTPOL):  ' 

READ  5,  INPUT3 

PRINT  *,'  Enter  name  of  output  file:  ' 

READ  5,  OUTPUT 
OPEN  (8, OUTPUT) 

WRITE  (8,19)  OUTPUT 

19  FORMAT ('Program  FRECl,  output  file  ',A) 

WRITE  (8,20)INPUT1 

IF(INDEP.EQ.2)  WRITE  (8,20)  INPUT2 
WRITE  (8,20)  INPUT3 

20  F0RHAT(' Input  file  - ',A) 

IF( INDEP. EQ.l)  WRITE  (8,1) 

WRITE  (8,1) 

WRITE  (8,21) 

21  FORMAT ( 'Rectangular  field  extends  from  -A  to  A In  x,  and  from  -B  t 
lo  B In  y. ' ) 

WRITE  (8,1) 

WRITE  (8,25) 

25  FORHAT('  TIN  - energy  of  Incident  proton  beam,  HeV') 

WRITE  (8,26) 

26  F0RHAT('  RGIN  - CSDA  range  at  energy  TIN,  g/cm2') 

WRITE  (8,27) 

27  FORHAT('  2A  - side  of  rectangular  field,  cm') 

WRITE  (8,28) 

28  FORHAT('  2B  - side  of  rectanglular  field,  cm') 

WRITE  (8,29) 

29  FORHAT('  X • x-dlstance  from  center  of  field,  cm') 

WRITE  (8,30) 

30  F0RMAT('  Y - Y-dlstance  from  center  of  field,  cm') 

WRITE  (8,31) 

31  FORMAT('’  LMAX  - number  of  (X,Y)  values') 

WRITE  (8,32) 

32  FORHAT('  NCASE  - number  of  depths') 

WRITE  (8,33) 

33  FORHAT('  REDUC  - reduction  factor') 

WRITE  (8,34) 


34 


35 


36 

37 

38 

39 


42 


43 


44 


45 

46 

47 

48 


100 

110 

115 

120 

125 


FORHAT('  DRAD  • absorbed  dose,  HeV/g') 

WRITE  (8,1) 

OPEN  (6, INPUTl) 

IF(INDEP.EQ.2)  OPEN  (7,INPUT2) 

READ  (6,5)  LINE 
READ  (6,5)  LINE 
READ  (6,«)  NCASEl 
READ  (6,«)  JHAXl 
DO  35  LN-1,8 
READ  (6,5)  LINE 
CONTINUE 

READ  (6,*)  TIN, RGIN 
READ  (6,5)  LINE 
NCASE2-0 
JHAX2-0 

IF( INDEP. EQ.l)  GO  TO  37 
READ  (7,5)  LINE 
READ  (7,5)  LINE 
READ  (7,*)  NCASE2 
READ  (7,*)  JMAX2 
DO  36  LN-1,5 
READ  (7,5)  LINE 
CONTINUE 

NCASE-NCASEl-)-NCASE2 
WRITE  (8,38) 

F0RHAT(9X,'TIN',8X,'RGIN',11X,'A',11X,'B','  LMAX  NCASE' 
WRITE  (8,39)  TIN, RGIN, A, B, LMAX, NCASE 
FORMAT(4F12.5,2I6) 

WRITE  (8,1) 

OPEN(ll, INPUTS) 

DO  42  LN-1,9 
READ  (11,5)  LINE 
CONTINUE 

READ  (11,*)  G1,G2,G3,G4,KHAX 
READ  (11,5)  LINE 
READ  (11,5)  LINE 
READ  (11,*)  (2D(K),K-1,KMAX) 

READ  (11,*)  LINE 

READ  (11,*)  (DEPOT(K),K=l,KMAX) 

CLOSE  (11) 

CALL  SC0FD( ZD , DEPOT , KMAX , AP , BP , CP , DP ) 

DO  250  NC-1, NCASE 

IF(NC-NCASE1)43,43,44 

JHAX-JHAXl 

READ  (6,5)  LINE 

READ  (6,*)  2R(NC),RHAX,SUM 

READ  (6,5)  LINE 

READ  (6,*)  (R(J),J-1,JMAX) 

READ  (6,5)  LINE 

READ  (6,*)  (F(J),J-1,JMAX) 

GO  TO  47 

JHAX-JHAX2 

READ  (7,5)  LINE 

READ  (7,*)  TIN,ZR(NC),RMAX 

READ  (7,5)  LINE 

READ  (7,*)  (R(J),J=1,JMAX) 

READ  (7,5)  LINE 

READ  (7,*)  (F(J),J-1,JMAX) 

JTAL-0 

DO  45  J-1,JHAX 

IF(F(J))46,46,45 

JTAL-JTAL+1 

JMAX-JTAL 

DO  48  J-1,JHAX 

FL(J)-LO6(2.0DO*PI*F(J)) 

CALL  SC0FD( R , FL , JHAX, ACOF , BCOF , CCOF , DCOF) 

DO  205  L-1,LHAX 
NHAX>NMX 
OLDVAL-2.0D0 
X-XP(L) 

Y-YP(L) 

IF(X.LE.A.AND.Y.LE.B)  GO  TO  110 
IF(X.LE.A.AND.Y.GT.B)  GO  TO  140 
IF(X.GT.A.AND.Y.LE.B)  GO  TO  160 
IF(X.GT.A.ANO.Y.GT.B)  GO  TO  180 
RHMAX=HIN(A-X,B-Y) 

IF(RMAX-RHMAX)115,115,120 
REDUC (L)-l.ODO 
GO  TO  200 

DRH-RHHAX/DBLE(NHAX-1 ) 

DO  125  N-1,NHAX 
RHO-DRH*DBLE(N-l) 

CALL  BSPOLD ( RHO , R , ACOF , BCOF , CCOF , DCOF , JMAX, RES ) 
GRAND(N)-RHO*EXP(RES) 

DELTA-DRH/3.0D0 

CALL  GRALD(DELTA,GRAND,NMAX, PARTI) 

RHHIN-HIN(A-X,B-Y) 

RHMAX-SQRT((X+A)**2+(Y+B)**2) 

RHHAX-HIN(RHHAX,RHAX) 

DRH-(RHMAX-RHMIN)/DBLE(NMAX-1) 

DO  130  N-1,NHAX 
RHO-RHHIN-t-ORH*DBLE(N-l) 

CALL  BSPOLD ( RHO , R , ACOF , BCOF , CCOF , DCOF , JHAX , RES ) 
FINT-RHO*EXP(RES) 

CALL  CORNER(RHO,A+X,B+Y,FRACl) 

CALL  CORN ER( RHO, A+X,B-Y,FRAC2) 
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CALL  C0RNER(RH0,A-X,B+Y,FRAC3) 

CALL  C0RNER(RH0,A-X,B-Y,FRAC4) 

130  GRAND(N)-FINT*(FRAC1+FRAC2+FRAC3+FRAC4) 
DELTA-DRH/3.000 

CALL  GRALD(DELTA, GRAND, NHM,PART2) 
REDUC(L)-PART1+PART2 
GO  TO  200 
140  RHHIN-Y-B 

I F ( RMAX-RHH IN ) 1 42 , 142 , 1 45 
142  REDUC(L}-0.0D0 
GO  TO  200 

145  RHHAX-SQRT((X+A)**2+(Y+B)**2) 

RHHAX-H IN (RHNAX , RHAX) 
DRH-(RHHAX-RHHIN)/DBLE(NHAX-1) 

DO  150  N-1,NHAX 
RHO-RHH IN+DRH*DBLE ( N- 1 ) 

CALL  BSPOLD (RHO , R , ACOF , BCOF , CCOF , DCOF , JMAX, RES ) 
FINT-RHO*EXP(RES) 

CALL  CORNER(RHO,A+X,Y+B,FRACl) 

CALL  CORNER(RHO,A-X,Y+B,FRAC2) 

CALL  CORNER(RHO,A+X,Y-=B,FRAC3) 

CALL  C0RNER(RH0,A-X,Y-B,FRAC4) 

150  GRAND(N)-FINT*(FRAC1+FRAC2-FRAC3-FRAC4) 
DELTA.DRH/3.0D0 

CALL  GRALD(OELTA, GRAND, NHAX,REDUC(L)) 

GO  TO  200 
160  RHHIN-X-A 

IF(RHAX-RHHIN)162,162,165 
162  REDUC(L)-O.ODO 
GO  TO  200 

165  RHHAX-SQRT((X+A)**2+(Y+B)**2) 
RHHAX-HIN(RHHAX,RHAX} 
DRH-(RHHAX-RHHIN)/DBLE(NHAX-1) 

DO  170  N-1,NHAX 
RHO-RHH IN+DRH*DBLE ( N- 1 ) 

CALL  BSPOLD (RHO , R , ACOF , BCOF , CCOF , DCOF , JHAX, RES ) 
FINT-RHO*EXP(RES) 

CALL  CORNER (RHO, X+A,B+Y,FRAC1) 

CALL  CORNER(RHO,X+A,B-Y,FRAC2) 

CALL  CORNER(RHO,X-A,B+Y,FRAC3) 

CALL  C0RNER(RH0,X-A,B-Y,FRAC4) 

170  GRAND(N)-FINT*(FRAC1+FRAC2-FRAC3-FRAC4) 
DELTA-DRH/3.0D0 

CALL  GRALD( DELTA, GRAND, NHAX , REDUC ( L ) ) 

GO  TO  200 

180  RHHIN-SQRT((X-A)**2+(Y-B)**2) 
IF(RHM-RHHIN)182,182,185 
182  REDUC (L)-O.ODO 
GO  TO  200 

185  RHHAX-SQRT((X+A)**2+(Y+B)**2) 
RHH/U(-HIN(RHHAX,RHAX} 
DRH-(RHHAX-RHHIN}/DBLE(NHAX-1} 

DO  190  N-1,NHAX 
RHO-RHH IN+DRH*DBLE ( N- 1 ) 

CALL  BSPOLD ( RHO , R , ACOF , BCOF , CCOF , DCOF , JHAX , RES ) 
FINT-RHO*EXP(RES) 

CALL  CORNER ( RHO, X+A,Y+B,FRAC1) 

CALL  CORNER ( RHO, X-A,Y+B,FRAC2) 

CALL  CORNER{RHO,X+A,Y-B,FRAC3) 

CALL  CORNER{RHO,X-A,Y-B,FRAC4) 

190  GRAND(N)-FINT*(FRAC1-FRAC2-FRAC3+FRAC4) 
DELTA-DRH/3.0D0 

CALL  GRALD( DELTA , GRAND , NHAX , REDUC ( L } ) 

200  DIFF-ABS(OLDVAL-REDUC(L)) 
IF(DIFF"EREDUCL)205,201,201 

201  OLDVAL-REDUC(L) 

IF(NHAX-2561)202,202,205 

202  NHAX-2*NHAX-1 
GO  TO  100 

205  CONTINUE 

CALL  BSPOLD ( ZR( NC ) , ZD , AP , BP , CP , DP , KHAX, DDDX (NC ) ) 
DO  206  L-1,LHAX 

206  DRAD(L)-DOOX(NC)«REDUC(L) 

WRITE  (8,210)  ZR(NC) 

210  FORHAT('Depth  as  fraction  of  RGIN  - ',F5.3) 

WRITE  (8,1) 

WRITE  (8,221) 

221  F0RHAT(11X,  'XMIX,  'Y'  ,7X,  'REDUC'  ,8X,  'DRAD' ) 

DO  223  L-1,LHAX 

WRITE  (8,222)  XP(L) ,YP(L) ,REDUC(L) ,DRAD(L) 

222  FORHAT(4F12.5) 

223  CONTINUE 
WRITE  (8,1) 

print  245,  ncase,nc 

245  format('  ncase  - ',i3,'  nc  - ',i3) 

250  CONTINUE 
STOP 
END 


6 A-HAX(SA,SB) 

B-HIN(SA,SB) 

IF(R-B)10, 10,20 
10  FRAC-0.25D0 
RETURN 

20  IF(R-A)30,30,40 
30  FRAC-0.25D0-ACOS(B/R)/TPI 
RETURN 

40  IF(R*R-A*A-B*B)50, 50,60 
50  FRAC-0 . 25D0- ( ACOS ( B/R ) +ACOS ( A/R ) ) /TPI 
RETURN 

60  FRAC-0. ODO 
RETURN 
END 


SUBROUTINE  CORNER(R,SA,SB,FRAC) 

27  Har  88.  Computes  rel .arciength  (fracton  of  2 pi)  for  point  in 
corner  of  rectangle  with  sides  SA  and  SB. 

IHPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

TPI-8.0D0*ATAN(1.0D0) 

IF(SA)60,60,5 
5 IF(SB)60,60,6 
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PROGRAM  FREC2 

15  May  93.  Written  by  Martin  J.  Berger,  NIST. 

Calculates  reduction  factors  and  absorbed-dose  values 
for  a rectangular  field.  FREC2  Is  similar  to  FRECl, 
but  carries  out  numerical  quadrature  with  respect 
to  X and  y. 

Center  of  coordinate  system  Is  at  x-0,  y»0. 

Rectangle  extends  from  -A  to  A In  x,  -B  to  B In  y. 

Results  are  calculated  a specified  set  of  (x,y)  values. 
Without  loss  of  generality,  It  Is  assumed  that  x and  y 
are  non-negative. 

The  external  functions  FUNl  and  FUN2  specify  the  upper  and  lower 
limits  of  y,  and  the  external  function  FUN3  specifies  a constant 
density  of  pencil  beams,  n(x,y)-l.  By  modifying  these  external 
functions,  the  FREC2  program  could  be  applied  to  fields  with 
arbitrary  shapes  and  specified  pencil-beam  density. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

EXTERNAL  FUN1,FUN2, FUN3 

DIMENSION  REDUC (401 ) , XP (401 ) , YP (401 ) , R(401 ) , F(401 ) , FL(401 ) , 

1 ACOF (401 ) , BCOF (401 ) , CCOF (401 ) , DC0F(401 ) , GRAND(501 ) , 

2 ZD(100),DEPOT(100),AP(100),BP(100),CP(100),DP(100),DRAD(401), 

3 ZR(401),DDDX(401), TERM (501) 

CHARACTER*16  INPUTl , INPUT2, INPUT3, OUTPUT, INFIL 
CHARACTER  LINE*80 
1 FORMAT (IH  ) 

5 FORMAT(A) 

PI-4.0D0*ATAN(1.0DO) 

PRINT  Rectangular  field  extends  from  -A  to  A In  x,  -B  to  B 1n 

lx' 

PRINT  Enter  values  of  A and  B (cm):  ' 

READ  *,  A,B 

PRINT  Options  for  specifying  (x,y )-va1ues: ' 

PRINT  1)  Enter  values  from  keyboard' 

PRINT  2)  Enter  values  from  file:' 

PRINT  *,'  Enter  choice:  ' 

READ  *,  INP 
IF(INP.EQ.2)  GO  TO  7 

PRINT  *,'  Enter  number  of  (x,y)  values:  ' 

READ  *,  LMAX 
DO  6 L>1,LMAX 

PRINT  Enter  x and  y (cm,  non-negative):  ' 

READ  *,  XP(L),YP(L) 

6 CONTINUE 
GO  TO  9 

7 PRINT  Enter  name  of  Input  file:  ' 

READ  5,  INFIL 

OPEN  (7, INFIL) 

READ  (7,*)  LMAX 

DO  8 L-1,LMAX 

READ  (7,*)  XP(L),YP(L) 

8 CONTINUE 
CLOSE (7) 

9 PRINT  Specification  of  Integration  grid:' 

PRINT  Enter  number  of  grid  points  In  x:  ' 

READ  *,  IGRD 

PRINT  Enter  number  of  grid  points  In  y:  ' 

READ  *,  JGRD 

PRINT  *,'  Enter  name  of  Input  file  1 (from  MORAD):  ' 

READ  5,  INPUTl 

PRINT  *,'  Options  for  depths:' 

PRINT  *,'  1)  No  depths  greater  than  98%  of  CSDA  range' 

PRINT  2)  Some  depths  greater  than  98%  of  CSDA  range' 

PRINT  *,'  Enter  choice:  ' 

READ  *,  INDEP 
IF(INDEP.EQ.l)  GO  TO  16 

PRINT  Enter  name  of  Input  file  2 (from  PTRAD):  ' 

READ  5,  INPUT2 

16  PRINT  *,'  Enter  name  of  Input  file  3 (from  PTPOL):  ' 

READ  5,  INPUT3 

PRINT  *,'  Enter  name  of  output  file:  ' 

READ  5,  OUTPUT 
OPEN  (8, OUTPUT) 

WRITE  (8,19)  OUTPUT 

19  FORMAT ('Program  FGEN,  output  file  ',A) 

WRITE  (8,20)INPUT1 
IF(INDEP.EQ.2)  WRITE  (8,20)  INPUT2 
WRITE  (8,20)  INPUT3 

20  FORMAT(' Input  file  - ',A) 

IF(INDEP.EQ.l)  WRITE  (8,1) 

WRITE  (8,1) 

WRITE  (8,21) 

21  FORMATf 'Rectangular  field  extends  from  -A  to  A In  x,  and  from  -B  t 
lo  B In  y.') 

WRITE  (8,1) 

WRITE  (8,25) 

25  FORMAT('  TIN  - energy  of  Incident  proton  beam,  MeV') 

WRITE  (8  26) 

26  FORMAT('  RGIN  » CSDA  range  at  energy  TIN,  g/cm2') 

WRITE  (8,27) 

27  FORMAT('  2A  - side  of  rectangular  field,  cm') 


■ side  of  rectanglular  field,  cm') 
number  of  grid  points  In  x') 
number  of  grid  points  In  y') 

■ x-dlstance  from  center  of  field,  cm') 

■ Y-dlstance  from  center  of  field,  cm') 

• number  of  (X,Y)  values') 

• number  of  depths' ) 

• number  of  grid  points  In  x for  numerical  quadrat 

■ nunber  of  grid  points  In  y for  numerical  quadrat 

• reduction  factor') 

■ absorbed  dose,  MeV/g') 


WRITE  (8,28) 

28  FORMAT ^ 2B 
WRITE  (8,29) 

29  FORMAT ('  IGRD 
WRITE  (8,30) 

30  FORMAT ('  JGRD 
WRITE  (8,31) 

31  FORMAT ( ' X 
WRITE  (8,32) 

32  FORMAT ('  Y 
WRITE  (8,33) 

33  FORMAT ('  LMAX 
WRITE  (8,34) 

34  FORMAT ^ NCASE 
WRITE  (8,35) 

35  FORMAT ('  IGRD 
lure' ) 

WRITE  (8,36) 

36  FORMAT ('  JGRD 
lure' ) 

WRITE  (8,37) 

37  FORMAT ^ REDUC 
WRITE  (8,38) 

38  FORMAT('  DRAD 
WRITE  (8,1) 

OPEN  (6, INPUTl) 

IF(INDEP.EQ.2)  OPEN  (7,INPUT2) 

READ  (6,5)  LINE 
READ  (6,5)  LINE 
READ  (6,*)  NCASEl 
READ  (6,*)  JMAXl 
DO  39  LN-1,8 
READ  (6,5)  LINE 

39  CONTINUE 

READ  (6,*)  TIN, RGIN 
READ  (6,5)  LINE 
NCASE2-0 
JMAX2-0 

IF( INDEP. EQ.l)  GO  TO  41 
READ  (7,5)  LINE 
READ  (7,5)  LINE 
READ  (7,*)  NCASE2 
READ  (7,*)  JMAX2 
DO  40  LN-1,5 
READ  (7,5)  LINE 

40  CONTINUE 

41  NCASE-NCASE1+NCASE2 
WRITE  (8,42) 

42  FORMAT(9X, 'TIN' ,8X, 'RGIN' ,11X, 'A' ,11X, 'B' ,2X, 'LMAX' ,1X, 'NCASE', 
1 2X, ' IGRD', 2X,' JGRD') 

WRITE  (8,43)  TIN, RGIN, A, B, LMAX, NCASE, IGRD, JGRD 

43  FORMAT(4F12.5,4I6) 

WRITE  (8,1) 

OPEN(ll,INPUT3) 

DO  44  LN-1,9 
READ  (11,5)  LINE 

44  CONTINUE 

READ  (11,*)  G1,G2,G3,G4,KMAX 
READ  (11,5)  LINE 
READ  (11,5)  LINE 
READ  (11,*)  (ZD(K),K-1,KMAX) 

READ  (11,*)  LINE 

READ  (11,*)  (DEP0T(K),K-1,KMAX) 

CLOSE  (11) 

CALL  SCOFD ( ZD , DEPOT , KMAX , AP , BP , CP , DP ) 

DO  250  NC-1, NCASE 
IF(NC-NCASE1)45,45,46 

45  JMAX-JMAXl 
READ  (6,5)  LINE 

READ  (6,*)  ZR(NC),RMAX,SUM 

READ  (6,5)  LINE 

READ  (6,*)  (R(J),J-1,JMAX) 

READ  (6,5)  LINE 

READ  (6,*)  (F(J),J-1,JMAX) 

GO  TO  49 

46  JMAX-JMAX2 
READ  (7,5)  LINE 

READ  (7,*)  TIN,ZR(NC),RMAX 

READ  (7,5)  LINE 

READ  (7,*)  (R(J),J«1,JMAX) 

READ  (7,5)  LINE 

READ  (7,*)  (F(J),J-1,JMAX) 

JTAL-0 

DO  47  J-1,JMAX 
IF(F(J))48,48,47 

47  JTAL-JTAL+1 

48  JMAX-JTAL 

49  DO  50  J-1,JMAX 

50  FL(J)-L0G(2.0D0*PI*F(J)) 

CALL  SCOFD( R, FL , JMAX, ACOF , BCOF , CCOF , DCOF ) 

DO  205  L-1,LMAX 
REDUC (L)-O.ODO 
XMIN— A 
XHAX*A 

DGX-(XMAX-XMIN)/DBLE(IGRD-1) 

DELTAX-DGX/3.0D0 
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DO  80  I-1,IGRD 
X-XMIN+DGX*DBLE(I-1) 
D6Y-(FUN2(X,B)-FUN1(X,B))/DBLE(JGRD-1) 
DELTAY-DGY/3.0D0 
DO  70  J-1,JGRD 
Y-FUN1(X,B)+DGY*DBLE(J-1) 
ARG-SQRT((XP(L)-X)**2+(YP(L)-Y)**2) 
IF(ARG-RHAX)60,60,55 
55  GRAND(J}-O.ODO 
GO  TO  70 

60  CALL  BSPOLD{ARG,R,ACOF,BCOF,CCOF,DCOF,JMAX,RES) 
GRAND(J)-FUN3(X,Y)*EXP{RES) 

70  CONTINUE 

CALL  GRALD( DELTAY , GRAND, JGRD, TERM ( I ) ) 

SO  CONTINUE 

CALL  GRALD(DELTAX,TERH,IGRD,SUHN) 

REDUC ( L ) -SUMH/ { 2 . ODO*PI ) 

205  CONTINUE 

CALL  BSPOLD(ZR(NC),ZD,AP,BP,CP,DP,KMAX,DODX(NC)) 
DO  206  L-1,LHAX 

206  ORAD(L}-DDDX(NC)«REDUC(L) 

WRITE  (8,210)  ZR(NC) 

210  FORHAT('Depth  as  fraction  of  RGIN  - ',F5.3) 

WRITE  (8,1) 

WRITE  (8,221) 

221  F0RHAT(11X,  'XMIX,  'Y'  ,7X,  'REDUC'  ,8X,  'DRAD' ) 

DO  223  L-1,LHAX 

WRITE  (8,222)  XP(L) ,YP(L) ,REDUC(L) ,DRAD(L) 

222  FORMAT (4F12. 5) 

223  CONTINUE 
WRITE  (8,1) 

print  245,  ncase,nc 

245  forniat('  ncase  - ',13,'  nc  - ',13) 

250  CONTINUE 
STOP 
END 

DOUBLE  PRECISION  FUNCTION  FUN1(X,B) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

FUNl— B 

RETURN 

END 

DOUBLE  PRECISION  FUNCTION  FUN2(X,B) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

FUN2-B 

RETURN 

END 

DOUBLE  PRECISION  FUNCTION  FUN3(X,Y) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-2) 

FUN3-1.0D0 

RETURN 

END 


PROGRAM  AXPLOT 


C 

C 24  July  93.  Written  by  Martin  J.  Berger,  NIST. 

C 

C Plots  reduction  factor  and  absorbed  dose  along  central  axis 
C of  circular  or  rectangular  field. 

C 

DIMENSION  Z ( 107) , D( 107) ,RED( 107 ) , A1 ( 107) , B1 { 107) , Cl ( 107 ) , D1  ( 107) , 
1 A2(107), B2{107),C2(107), 02(107), X{801),Y1(801),Y2(801) 

CHARACTER  INPUT»16,LINE*80,TITLE*80 
DATA  JHAX/801/ 

5 FORMAT(A) 

PRINT  Options  for  Input:' 

PRINT  1)  Circular  field' 

PRINT  2)  Rectangular  field' 

PRINT  *,'  Enter  choice:  ' 

READ  *,  IN 

PRINT  Enter  name  of  input  file:  ' 

READ  5,  INPUT 

PRINT  ♦,'  Enter  title  of  plot:  ' 

READ  5,  TITLE 
OPEN  (7, INPUT) 

GO  TO  (10,40),  IN 
10  DO  20  LN-1,12 
READ  (7,5)  LINE 
20  CONTINUE 

READ  (7,«)  G1,G2,G3,IG,NHAX 
DO  30  LN-1,4 
READ  (7,5)  LINE 
30  CONTINUE 

READ  (7,*)  (Z(N),N-1,NMAX) 

READ  (7,5)  LINE 

READ  (7,*)  (RE0(N),N-1,NMAX) 

READ  (7,5)  LINE 

READ  (7,*)  (D(N),N-1,NMAX) 

CLOSE  (7) 

GO  TO  70 

40  DO  50  LN-1,18 
READ  (7,5)  LINE 
50  CONTINUE 

READ  (7,*)  G1,G2,G3,G4,G5,IG,NNAX 
DO  60  LN-1,4 
READ  (7,5)  LINE 
60  CONTINUE 

READ  (7,*)  (Z(N),N-1,NMAX) 

READ  (7,5)  LINE 

READ  (7,*)  (RED(N),N-1,NNAX) 

READ  (7,5)  LINE 

READ  (7,«)  (0(N),N-1,NMAX) 

CLOSE  (7) 

70  CALL  SC0F(Z,RED,NNAX,A1,B1,C1,D1) 

CALL  SCOF(Z,D,NNAX,A2,B2,C2,D2) 

0X-Z(NMAX)/REAL(JHAX-1) 

DO  80  J-1,JHAX 
X(J)-0X*REAL(J-1) 

CALL  BSP0L(X(J),Z,A1,B1,C1,D1,NNAX,Y1(J)) 

CALL  BSPOL(X(J),Z,A2,B2,C2,D2,NHAX,Y2(J)) 

80  CONTINUE 

CALL  SETDV('HPG') 

100  CALL  LINLOG(0,0,1) 

CALL  SIDTEX(TITLE,l,'2/r_o',l, 'Reduction  Factor', 1,'  ',!) 

CALL  PLAC(1. 0,0. 5,0. 0,1.0) 

CALL  SETLIH(2,3,0.0,1.0) 

CALL  HOWPLT(0,1,1) 

CALL  CURV(JMAX,X,Y1} 

CALL  VG 

CALL  SIDTEX('  ', 1, 'z/r_o' ,1, 'Absorbed  Dose,  MeV/g' ,1, ' ',!) 

CALL  PLAC(0.5,0.0,0.0,1.0) 

CALL  H0WPLT(0,1,1) 

CALL  CURV(JHAX,X,Y2) 

CALL  VG 

IF(LOOPIN().Eq.l)  GO  TO  100 

STOP 

END 
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PROGRAM  RADPLOT 


C 

C 24  July  93.  Written  by  Martin  J.  Berger,  NIST. 

C 

C Plots  reduction  factor  and  absorbed  dose  at  a fixed  depth, 

C as  functions  of  the  radial  distance  from  the  center  of  the  field. 

DIMENSION  RHO(201), RED(201), 0(201), A1(201),B1(201),C1(201), 

1 Dl(201), A2(201),B2(201),C2(201), 02(201), X(801),Y1(801), 

2 Y2(801) 

CHARACTER  INPUT*16,LINE*80,TITLE*80 
DATA  JMAX/801/ 

5 FORHAT(A) 

PRINT  Options  for  Input:' 

PRINT  1)  Circular  field' 

PRINT  *,'  2)  Rectangular  field' 

PRINT  Enter  choice:  ' 

READ  IN 

PRINT  *,'  Enter  nane  of  Input  file:  ' 

READ  S,  INPUT 

PRINT  *,'  Enter  title  of  plot:  ' 

READ  5,  TITLE 
OPEN  (7, INPUT) 

GO  TO  (10,40),  IN 
10  DO  20  LN-1,12 
READ  (7,5)  LINE 
20  CONTINUE 

READ  (7,*)  G1,G2,63,LHAX 
READ  (7,5) 

READ  (7,5) 

READ  (7,*)  (RH0(L),L-1,LMAX) 

READ  (7,5)  LINE 
READ  (7,5)  LINE 
READ  (7,*)  (RED(L),L-1,LHAX) 

READ  (7,5)  LINE 

READ  (7,*)  (D(L),L-1,LMAX) 

GO  TO  60 

40  DO  50  LN-1,18 
READ  (7,5)  LINE 
50  CONTINUE 

READ  (7,«)  G1,G2,G3,G4,G5,LNAX 

READ  (7,5)  LINE 

READ  (7,5)  LINE 

READ  (7,*)  (RHO(L),L-l,LHAX) 

READ  (7,5)  LINE 
READ  (7,5)  LINE 
READ  (7,*)  (RED(L),L-1,LHAX) 

READ  (7,5)  LINE 

READ  (7,*)  (D(L),L-1,LHAX) 

60  CLOSE  (7) 

CALL  SC0F(RH0,RED,LHAX,A1,B1,C1,D1) 

CALL  SCOF(RHO,D,LHAX,A2,B2,C2,D2) 

DX-RHO( LMAX ) /REAL ( JMAX- 1 ) 

DO  70  J-1,JHAX 
X(J}-DX«REAL(J-1) 

CALL  BSP0L(X(J),RH0,A1,B1,C1,01,LNAX,Y1(J)) 

CALL  BSPOL(X(J),RHO,A2,B2,C2,02,LHAX,Y2(J)) 

X(J)-10.0*X(J) 

70  CONTINUE 

CALL  SET0V('HP6') 

100  CALL  LINLOG(0,0,1) 

CALL  SIDTEX(TITLE,l,'r, mm', 1, 'Reduction  Factor', 1,'  ',!) 

CALL  PLAC(1.0,0.5,0.0,1.0) 

CALL  SETLIH(2,3,0.0,1.0) 

CALL  HOWPLT(0,1,1) 

CALL  CURV(JHAX,X,Y1) 

CALL  VG 

CALL  SIDTEX('  ',1,'r,  nn' ,1, 'Absorbed  Dose,  HeV/g',1,'  ',!) 

CALL  PLAC(0. 5, 0.0,0. 0,1.0) 

CALL  H0WPLT(0,1,1) 

CALL  CURV(JMAX,X,Y2} 

CALL  VG 

IF(LOOPIN().EQ.l)  GO  TO  100 

STOP 

END 
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CALF-COS(ALF) 

SALF-SIN(ALF) 

DO  20S  L-l.LHAX 
NHAX-NHX 
OLOVAL-Z.ODO 
X-CALF*DS*DBLE(L-1) 

Y-SALF*DS*DBLE(L-1) 

100  IF(X,LE.A.AMO.Y.LE.B)  GO  TO  110 
IF(X.LE.A.ANO.Y.GT.B)  GO  TO  140 
IF{X.GT.A.AND.Y.tE.B)  GO  TO  160 
IF(X.GT.A.ANO.Y.GT.B)  GO  TO  180 
110  RHMAX-HIN(A-X,B-Y) 

IF(RMAX-RHHAX)115,115,120 
115  RE0UC(L)-1.0D0 
GO  TO  200 

120  ORH-RHHAX/OBLE(NHAX-l) 

DO  125  N-1,NHAX 
RHO-ORH*DBLE(N-l) 

CALL  BSPOLD(RHO,R,ACOF,BCOF,CCOF,DCOF,JHAX,RES) 
125  6RAND(N)-RHO«EXP(RES) 

0ELTA-0RH/3.0D0 

CALL  GRALD( DELTA, GRAND, HHAX, PARTI) 
RHHIN-MIN(A-X,B-Y) 

RHMAX-SQRT({X+A)**2+(Y+B)**2) 

RHNAX-N IN ( RHMAX . RHAX ) 
DRH-(RHMAX-RHHIN)/DBLE(NMAX-1 ) 

DO  130  N-1,NHAX 
RH0-RHNIN4^0RH*DBLE(N-1  ) 

CALL  BSPOLD ( RHO , R , ACOF , BCOF , CCOF , OCOF , JMAX, RES ) 
FINT-RHO*EXP(RES) 

CALL  C0RNER(RH0,A4^X,B-^Y,FRAC1) 

CALL  C0RNER(RH0,A+X,B-Y,FRAC2) 

CALL  C0RNER(RH0,A-X,B+Y,FRAC3) 

CALL  C0RNER(RH0,A-X,B-Y,FRAC4) 

130  6RAND(N)-FINT*(FRAC1+FRAC2+FRAC3+FRAC4) 
OELTA-ORH/3.000 

CALL  GRALD( DELTA, GRAND. NHAX.PART2) 
REDUC(L)-PART1+PART2 
GO  TO  200 
140  RHHIN-Y-B 

I F ( RMAX-RHM IN ) 142 , 142 , 1 45 
142  REOUC(L)-O.ODO 
GO  TO  200 

145  RHHAX-S(}RT((X+A)**2+(Y+B)**2) 
Rra<AX-NIN(RHHAX,RHAX) 
DRH-(RHMAX-RHNIN)/DBLE(NNAX-1 ) 

DO  150  N-1,NHAX 
RHO-RHNIN^ORH*OBLE(N-l) 

CALL  BSPOLO( RHO, R.ACOF, BCOF, CCOF, DCOF,JHAX,RES) 
FINT-RHO* EXPIRES) 

CALL  C0RNER(RH0,A>X,Y4^B,FRAC1) 

CALL  C0RNER(RH0,A-X,Y4^B,FRAC2) 

CALL  C0RNER(RH0,Aa,Y-B,FRAC3) 

CALL  C0RNER(RH0,A-X,Y-B,FRAC4) 

150  GRAN0(H)-FINT*(FRAC1+FRAC2-FRAC3-FRAC4) 
0ELTA-ORH/3.0D0 

CALL  GRALD(DELTA,GRANO,NHAX,REOUC(L)) 

GO  TO  200 
160  RHHIN-X-A 

IF(RNAX-RHHIN)162,162,165 
162  REDUC(L) -0.000 
GO  TO  200 

165  RHMAX-SQRT((X+A)««2+(Y+B)**2) 
RHHAX-NINIRHHAX.RHAX) 
DRH-<RHHAX-RHHIN)/DBLE(NHAX-1 ) 

DO  170  N-1,NHAX 
RHO-RHN IN40RH*DBLE ( N - 1 ) 

CALL  BSPOLO(RHO,R, ACOF, BCOF, CCOF, OCOF, JHAX,RES) 
FINT-RHO«EXP(RES) 

CALL  CORNER(RHO,X«A,B^Y,FRACl) 

CALL  C0RNER(RH0,XM,B-Y,FRAC2) 

CALL  C0RNER(RH0,X-A,B^Y,FRAC3) 

CALL  C0RNER(RH0,X-A,B-Y,FRAC4) 

170  GRAN0(N)-FINT«(FRAC1+FRAC2-FRAC3-FRAC4) 
DELTA-ORH/3.000 

CALL  GRALO(OELTA,GRAND,NHAX,REDUC(L)) 

GO  TO  200 

180  RHMIN-SqRT{(X-A)**2+(Y-B)**2) 

I F( RMAX-RHNIN ) 182 , 182 , 185 
182  REOUC(L) -0.000 
GO  TO  200 

185  RHHAX-SQRTI(X+A)**2+(Y+B)**2) 

RHNAX-H IN ( RHHAX , RHAX ) 

0RH-( RHHAX-RHNIN ) /DBLE( NHAX-1 ) 

DO  190  N-1,NNAX 
RHO-RHN  IN4^DRH*DBLE(  N- 1 ) 

CALL  BSPOLOIRHO, R.ACOF, BCOF, CCOF, OCOF, JHAX.RES) 
FINT-RHO* EXPIRES) 

CALL  C0RNERIRH0,X^A,Y4B,FRAC1) 

CALL  CORNER I RHO, X-A,Y48,FRAC2) 

CALL  C0RNERIRH0,X^A,Y-B,FRAC3) 

CALL  C0RNER(RH0,X-A,Y-B,FRAC4) 

190  GRA»IN)-FINT*IFRACl-FRAC2-FRAC3+FRAC4) 
DELTA-ORH/3.000 

CALL  GRALOI DELTA, GRAND. NMAX.REOUCID) 

200  OIFF-ABSIOLOVAL-REDUC(L)) 


IFIDIFF-EREOUCL)205,201,201 

201  OLOVAL-REDUCIL) 

I F I NHAX-256 1 ) 202 , 202 . 205 

202  NMAX-2*NHAX-1 
GO  TO  100 

205  CONTINUE 

CALL  BSPOLOIZRINC),ZO,AP,BP,CP,DP,KHAX,OOOXINC)) 

DO  206  L-l.LHAX 

206  DRAOIL)-OOOXINC)«REDUCIL) 

IFILHAX.NE.l)  WRITE  18,210)  ZR(NC) 

210  FORHATI'Oepth  as  fraction  of  R6IN  - ',F5.3) 

IF(LHAX.EQ.l)  REDINO-REDUCIl) 

IFILHAX.EQ.l)  0AXINC)-DRA0I1) 

IFILHAX.NE.l)  WRITE  18,220) 

220  FORHATI 'Reduction  factors') 

IFILHAX.NE.l)  WRITE  18,230)  (REDUCIL) ,L-1,LHAX) 

230  F0RHATI8F12.6) 

IFILHAX.NE.l)  WRITE  18,240) 

240  FORHATI 'Absorbed  dose,  HeV/g,  from  a bean  of  1 proton/cn2') 
IF(LHAX.NE.l)  WRITE  18,230)  (DRAOIL), L-l.LHAX) 
print  245,  ncase.nc 

245  fomttr  ncase  - ',13,'  nc  - ',13) 

250  CONTINUE 

IF(LHAX.NE.l)  GO  TO  290 
WRITE  18,260)  SU) 

260  FORHATI 'Radial  distance  (cm)  from  center  of  field  - ',F5.3) 
WRITE  (8,1) 

WRITE  (8,270) 

270  FORHATI 'Depths,  In  units  of  RGIN') 

WRITE  (8,280)  (ZR(NC) ,NC-1, NCASE) 

280  F0RHATI8F12.6) 

WRITE  (8,220) 

WRITE  (8,280)  (RED(NC) ,NC-1, NCASE) 

WRITE  (8,240) 

WRITE  (8,280)  (DAX(NC),NC-1, NCASE) 

290  STOP 
END 

SUBROUTINE  CORNERIR,SA,SB,FRAC) 

C 27  Har  88.  Computes  rel .arc length  (fracton  of  2 pi)  for  point 
C comer  of  rectangle  with  sides  SA  and  SB. 

IHPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

TPI-8.000*ATAN( 1.000) 

IF(SA)60,60,S 

5 IF(S8)60,60,6 

6 A-HAX(SA,SB) 

B-NINISA,SB) 

IF(R-B)10,10,20 

10  FRAC-0.2500 
RETURN 

20  IF(R-A)30,30,40 
30  FRAC-0.25D0-ACOS(B/R)/TPI 
RETURN 

40  IF(R*R-A*A-8*B)50,50,60 
50  FRAC-0.25DO-(ACOS(B/R)+ACOS(A/R))/TPI 
RETURN 

60  FRAC-O.ODO 
RETURN 
END 
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PROGRAM  FREC 

15  May  93.  Written  by  Martin  J.  Berger,  NIST. 

Calculates  reduction  factors  and  absorbed-dose  values 
for  a rectangular  field. 

Center  of  coordinate  system  Is  at  x-0,  y-0. 

Rectangle  extends  from  -A  to  A In  x,  -B  to  B In  y. 

Results  are  calculated  as  functions  of  the  distance 
from  the  center  of  the  field,  along  a line  that  makes 
an  angle  alpha  (between  0 and  90  degrees)  with  respect 
to  the  x-axIs. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  REDUC(401},S(401),R(401),F(401),FL(401), 

1 ACOF ( 401 ) , BCOF ( 401 ) , CCOF ( 401 ) , OCOF ( 401 ) , GRAND( 5121 ) , 

2 ZD(100),OEPOT(100),AP(100),BP(100),CP(100),DP(100),DRAD(401), 

3 ZR(4Ol),DODX(401),RE0(4Ol),0AX(401) 

CHARACTER«16  INPUTl , INPUT2, INPUTS , OUTPUT, INFIL 
CHARACTER  LINE*80 

DATA  NMX/161/,EREDUCL/1.0D-04/ 

I FORMAT (IH  ) 

5 FORMAT(A) 

PI-4.000*ATAN(1.000) 

RAOFAC-PI/180.0D0 

PRINT  Rectangular  field  extends  from  -A  to  A In  x,  -B  to  B In 
ly' 

PRINT  Enter  values  of  A and  B (cm):  ' 

READ  *,  A,B 

PRINT  *,  ' Line  goes  through  a corner  of  the  rectangle  (1-yes,  2-n 
lo):  ' 

READ  *,  JVEC 
IF(JVEC.HE.2)  GO  TO  10 

PRINT  Enter  angle  alpha  with  respect  to  x-axIs  (degrees):  ' 
READ  *,  ALPHA 
ALF-ALPHA*RAOFAC 
GO  TO  11 
10  ALF-ATAN{B/A) 

ALPHA-ALF/RAOFAC 

II  PRINT  Options  for  specifying  radial  distances:' 

PRINT  1)  Get  set  of  values  from  file' 

PRINT  2)  Get  set  of  values  from  keyboard' 

PRINT  3)  Specify  set  of  values  In  terms  of  maximum  distance' 

PRINT  *,'  and  number  of  distances,  LHAX:  ' 

PRINT  *,'  Enter  choice:  ' 

READ  *,  INS 

GO  TO  (12,13,14),  INS 

12  PRINT  Enter  name  of  file  with  distance  values:  ' 

READ  S,  INFIL 

OPEN  (7. INFIL) 

READ  (7,*)  LMAX 

READ  (7,*)  (S(L),L-1,LHAX) 

CLOSE  (7) 

GO  TO  IS 

13  PRINT  *,'  Enter  number  of  radial  distances  (no  greater  than  401): 
1' 

READ  *,  LMAX 

PRINT  Enter  distances  (cm):  ' 

READ  ♦,  (S(L),L-1,LMAX) 

GO  TO  15 

14  PRINT  Enter  maximum  radial  distance  (cm):  ' 

READ  *,  SMAX 

LMAX-l 

IF(SHAX.LE. 0.000)  GO  TO  15 

PRINT  Enter  nunber  of  radial  distances  (no  greater  than  401): 
1' 

READ  *,  LHAX 

15  PRINT  *,'  Enter  name  of  Input  file  1 (from  HORAD):  ' 

READ  5,  INPUTl 

PRINT  *,'  Options  for  depths:' 

PRINT  *,'  1)  No  depths  greater  than  98S  of  CSOA  range' 

PRINT  *,'  2)  Some  depths  greater  than  98X  of  CSOA  range' 

PRINT  •,'  Enter  choice:  ' 

READ  *,  INDEP 
IF(INOEP.EQ.l)  GO  TO  16 

PRINT  *,'  Enter  name  of  Input  file  2 (from  PTRAD):  ' 

READ  5,  INPUT2 

16  PRINT  ♦,'  Enter  name  of  Input  file  3 (from  PTPOL):  ' 

READ  5,  INPUTS 

PRINT  Enter  name  of  output  file:  ' 

READ  5,  OUTPUT 
IF(INS.LT.3)  GO  TO  18 
S(l)-0.0 

IF(LHAX.LE.l)  GO  TO  18 
0S-SMAX/0BLE( LMAX-l) 

DO  17  L-2,LHAX 

17  S(L)-0S*0BLE(L-1) 

18  OPEN  (8,OUTPUT) 

WRITE  (8,19)  OUTPUT 

19  format! 'Program  FREC,  output  file  ' ,A) 

WRITE  (8, 20) INPUTl 
IF(INDEP.EQ.2)  WRITE  (8,20)  INPUT2 
WRITE  (8,20)  INPUTS 

20  F0RHAT( ' Input  file  - ',A) 


IF(INDEP.EO.l)  WRITE  (8,1) 

WRITE  (8,1) 

WRITE  (8,25) 

25  format!'  TIM  - energy  of  Incident  proton  beam,  MeV') 

WRITE  (8,26) 

26  format!'  RGIN  - CSOA  range  at  energy  TIN,  g/cm2') 

WRITE  (8,27) 

27  format!'  2A  - side  of  rectangular  field,  cm') 

WRITE  (8,28) 

28  format!'  2B  - side  of  rectanglular  field,  cm') 

WRITE  (8f29) 

29  F0RHAT('’  ALPHA  - angle  (deg)  defined  below') 

WRITE  (8,30) 

30  format!'  LHAX  - number  of  radial  distances  from  center  of  field' 
1) 

WRITE  (8,31) 

31  format!'  NCASE  - number  of  depths') 

WRITE  (8,1) 

WRITE  (8,32) 

32  format! 'Rectangular  field  extends  from  -A  to  A In  x,  and  from  -B  t 
lo  B In  y. ' ) 

WRITE  (8,33) 

33  format! 'Line  along  which  reduction  factor  Is  calculated  starts  at 
Ithe  origin' ) 

WRITE  (8,34) 

34  format! 'and  makes  angle  ALPHA  (deg)  with  respect  to  x-ax1s.') 

WRITE  (8,1) 

OPEN  (6, INPUTl) 

IF(INDEP.EQ.2)  OPEN  (7,INPUT2) 

READ  (6,5)  LINE 
READ  (6,5)  LINE 
READ  (6,«)  NCASEl 
READ  (6,*)  JMAXl 
DO  35  LN-1,8 
READ  (6,5)  LINE 

35  CONTINUE 

READ  (6,*)  TIM,RGIN 
READ  (6,5)  LINE 
NCASE2-0 
JHAX2-0 

IF( INDEP. EP.l)  GO  TO  37 
READ  (7,5)  LINE 
READ  (7,5)  LINE 
READ  (7,*)  NCASE2 
READ  (7,*)  JHAX2 
DO  36  LN-1,5 
READ  (7,5)  LINE 

36  CONTINUE 

37  NCASE-NCASE14NCASE2 
WRITE  (8,38) 

38  F0RHAT(9X,'TIN',8X,'RGIN',11X,'A',11X,'B',7X,'ALPHA','  LMAX  NCASE 
1') 

WRITE  (8,39)  TIN,RGIN,A,B,ALPHA,LHAX,NCASE 

39  F0RMAT(5F12.5,2I6) 

IF(LMAX.NE.l)  WRITE  (8,1) 

IF(LMAX.NE.l)  WRITE  (8,40) 

40  F0RMAT( 'Radial  distances  from  center  of  field  (cm)') 

IF(LMAX.NE.l)  WRITE  (8,41)  (S(L),L-1,LMAX) 

41  FORMAT (6F12. 5) 

0PEN( 11, INPUTS) 

DO  42  LN-1,9 
READ  (11,5)  LINE 

42  CONTINUE 

READ  (11,*)  G1,62,G3,G4,KMAX 
READ  (11,5)  LINE 
READ  (11,5)  LINE 
READ  (11,*)  (ZD(K),K-1,KMAX) 

READ  (11,*)  LINE 

READ  (11,*)  (DEPOT(K),K-l,KMAX) 

CLOSE  (11) 

CALL  SCOFD(ZD,OEPOT,KMAX,AP,BP,CP,DP) 

DO  250  NC-1, NCASE 
IF(NC-NCASE1)43, 43,44 

43  JNAX-JMAXl 
READ  (6,5)  LINE 

READ  (6,*)  ZR(NC),RNAX,SUH 

READ  (6,5)  LINE 

READ  (6,*)  (R(J),J-1,JHAX) 

READ  (6,5)  LINE 

READ  (6,*)  (F(J),J-1,JHAX) 

GO  TO  47 

44  JMAX-JHAX2 
READ  (7,5)  LINE 

READ  (7,«)  TIN,ZR(NC),RHAX 

READ  (7,5)  LINE 

READ  (7,*)  (R(J),J-1,JHAX) 

READ  (7,5)  LINE 

READ  (7,*)  (F(J),J-1,JHAX) 

JTAL-0 

DO  45  J-1,JHAX 
IF(F(J))46,46,45 

45  JTAL-JTAL^l 

46  JMAX-JTAL 

47  DO  48  J-1,JHAX 

48  FL(J)-L0G(2.0D0*PI*F(J)) 

CALL  SC0F0( R, FL , JHAX, ACOF , BCOF , CCOF , OCOF ) 
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PROGRAM  AXPLOT 


C 

C 24  Juljr  93.  Written  by  Martin  J.  Berger,  NIST. 

C 

C Plots  reduction  factor  and  absorbed  dose  along  central  axis 
C of  circular  or  rectangular  field. 

C 

DIMENSION  Z(1071,D(107),RED(107),A1(107),B1(107),C1(107),D1(107K 
1 A2(107), B2(107),C2(107), 02(107), X(801),Y1(801),Y2{801) 

CHARACTER  INPUT*16, LINE*80,TITLE*80 
DATA  JHAX/801/ 

5 FORMAT! A) 

PRINT  Options  for  Input:' 

PRINT  1)  Circular  field' 

PRINT  ♦,'  2)  Rectangular  field' 

PRINT  Enter  choice:  ' 

READ  *,  IN 

PRINT  Enter  name  of  Input  file:  ' 

READ  5,  INPUT 

PRINT  Enter  title  of  plot:  ' 

READ  5,  TITLE 
OPEN  (7, INPUT) 

GO  TO  (10,40),  IN 
10  DO  20  LN-1,12 
READ  (7,5)  LINE 
20  CONTINUE 

READ  (7,*)  G1,G2,G3,IG,NMAX 
DO  30  LN-1,4 
READ  (7,5)  LIME 
30  CONTINUE 

READ  (7,*)  (Z(N),N-1,NHAX) 

READ  (7,5)  LINE 

READ  (7,*)  (RED(N),N-1,NHAX) 

READ  (7,5)  LINE 

READ  (7,*)  (0(N),N-1,NMAX) 

CLOSE  (7) 

GO  TO  70 

40  DO  50  LN-1,18 
READ  (7,5)  LINE 
50  CONTINUE 

READ  (7,«)  G1,G2,G3,G4,G5,IG,NMAX 
DO  60  LN-1,4 
READ  (7,5)  LINE 
60  CONTINUE 

READ  (7,«)  (Z(N),N-1,NHAX) 

READ  (7,5)  LINE 

READ  (7,*)  (RED(N),N-1,NHAX) 

READ  (7,5)  LINE 

READ  (7,«)  (D(N),N-1,NHAX) 

CLOSE  (7) 

70  CALL  SC0F(Z,RE0,NHAX,A1,B1,C1,01) 

CALL  SC0F(Z,0,NMAX,A2,B2,C2,02) 

DX-Z(NHAX)/REAL(JHAX-1) 

DO  80  J-1,JHAX 
X(J)-0X*REAL(J-1) 

CALL  BSP0L(X(J),Z,A1,B1,C1,01,NHAX.Y1(J)) 

CALL  BSPOL(X(J),2,A2,82,C2,D2,NMAX,Y2(J)) 

80  CONTINUE 

CALL  SETOV('HPG') 

100  CALL  L1NLOG(0.0,1) 

CALL  SIDTEX(TITLE,l,'z/r  o' , 1, 'Reduction  Factor', 1,'  ',1) 

CALL  PLAC(1.0,0.5,0.0,1.0) 

CALL  SETLIH(2,3,0.0,1.0) 

CALL  HOWPLT(0,1,1) 

CALL  CURV(JHAX.X,Y1) 

CALL  VG 

CALL  SIDTEX('  ', 1, '2/r_o' ,1, 'Absorbed  Dose,  MeV/g' , 1, ' ',!) 

CALL  PLAC(0.5,0.0,0.0,1.0) 

CALL  HOWPLT(0,1,1) 

CALL  CURV(JHAX,X,Y2) 

CALL  VG 

IF(LOOPIN().Eq.l)  GO  TO  100 

STOP 

END 
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PROGRAM  RAOPLOT 


C 

C 24  July  93.  Written  by  Martin  J.  Berger,  NIST. 

C 

C Plots  reduction  factor  and  absorbed  dose  at  a fixed  depth, 

C as  functions  of  the  radial  distance  from  the  center  of  the  field. 

DIMENSION  RHO(201),REO(201),D(201),A1(201),B1(201),C1(201), 

1 Dl(201), A2(201),B2(201),C2(201), 02(201), X{801),Y1(801), 

2 Y2(801) 

CHARACTER  INPUT*16, LINE*80,TITLE*80 
DATA  JHAX/801/ 

5 FORHAT(A} 

PRINT  Options  for  Input:' 

PRINT  1)  Circular  field' 

PRINT  2)  Rectangular  field' 

PRINT  Enter  choice:  ' 

READ  «,  IN 

PRINT  *,'  Enter  name  of  Input  file:  ' 

READ  S,  INPUT 

PRINT  Enter  title  of  plot:  ' 

READ  5,  TITLE 
OPEN  (7, INPUT) 

GO  TO  (10,40),  IN 
10  DO  20  LN-1,12 
READ  (7,5)  LINE 
20  CONTINUE 

READ  (7,*)  G1,G2,G3,LHAX 
READ  (7,5) 

READ  (7,5) 

READ  (7,*)  (RHO(L),L-l,LHAX) 

READ  (7,5)  LINE 
READ  (7,5)  LINE 
READ  (7,*)  (RED(L),L-1,LHAX) 

READ  (7,5)  LINE 

READ  (7,*)  (D(L),L-1,LHAX) 

GO  TO  60 

40  DO  50  LN-1,18 
READ  (7,5)  LINE 
50  CONTINUE 

READ  (7,*)  61,G2,G3,G4,6S,LMAX 

READ  (7,5)  LINE 

READ  (7,5)  LINE 

READ  (7,*)  (RHO(L),L-l,LHAX) 

READ  (7,5)  LINE 
READ  (7,5)  LINE 
READ  (7,*)  (RED(L),L-1,LHAX) 

READ  (7,5)  LINE 

READ  (7,*)  (D(L),L-1,LMAX) 

60  CLOSE  (7) 

CALL  SC0F(RH0,RED,LHAX,A1,B1,C1,D1) 

CALL  SC0F(RH0,D,LNAX,A2,B2,C2,02) 

DX-RHO( LNAX ) /REAL ( JHAX- 1 ) 

DO  70  J-1,JHAX 
X(J)-DX*REAL(J-1) 

CALL  BSPOL(X(J),RHO,A1,B1,C1,01,LMAX,Y1(J)) 

CALL  BSP0L(X(J),RH0,A2,B2,C2,02,LNAX,Y2(J)) 

X(J)-10.0*X(J) 

70  CONTINUE 

CALL  SETDV('HPG') 

100  CALL  LINLOG(0,0,1) 

CALL  SIDTEX(TITLE,l,'r, mm', 1, 'Reduction  Factor', 1,'  '.1) 

CALL  PLAC(1. 0,0. 5,0. 0,1.0) 

CAU  SETLIH(2,3,0.0,1.0) 

CALL  HOWPLT(0,1,1) 

CALL  CURV(JHAX,X,Y1) 

CALL  VG 

CALL  SIDTEX('  1, 'r,  mm' ,1, 'Absorbed  Dose,  MeV/g' , 1, ' ',!) 

CALL  PLAC(0. 5, 0.0,0. 0,1.0} 

CALL  H0WPLT(0,1,1} 

CALL  CURV(JMAX,X,Y2} 

CALL  VG 

IF(LOOPIN().EQ.l)  GO  TO  100 

STOP 

END 
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